//Solve bioreactor modelling problem Ahmad et al 2011.
dirbioreactor = get_absolute_file_path("bioreactor-PRODA.sce")
function F = ff(t, y)//Function to solve the differential equation
mumax = 0.084;// per hour
ks = 213.6 ;//g/l
yxs = 0.136;//g/g
ypx = 4.913;//g/g
yps = 0.6682;//g/g
qp = 0.4277;// per hour

    F = zeros(1,3);//initialize
    mu = (mumax * y(2))/(ks + y(2));//equation (1)
    F(1) = mu * y(1);//equation (2)
    F(2) = - mu * y(1)/yxs;//equation (3)
    F(3) = qp * y(1);//equation (4)
endfunction
y0 = [2.12  240     0];//initial value of y
t0 = 0;//initial time
t = [0:0.01:80];//values of t for which the functions are to be evaluated
y=ode(y0,t0,t,ff);//solve the differential equation
k = length(y);//number of elements of y
//y1 are the odd elements, while y2 are the even number elements
y1i = 1:3:k;//select index of  y1
y1 = y(y1i);
y2i = 2:3:k;//select index of y2 starting at 2
y2 = y(y2i);
y3i = 3:3:k
y3 = y(y3i);//select index of y3 starting at 3
result = [t' y1' y2' y3'];//arrange in a table
//disp(result);
scf(0);//create new graphic window
plot(t, y1);//plot of biomass concentration with time
xtitle("Plot of biomass concentration", "Time, hours", "Biomass concentration, g/g");
myfile = "bioreactor-biomass.png";//file name
filename = dirbioreactor + myfile;
xs2png(0, filename);//save file as PNG
scf(1);//create new graphic window
plot(t, y2);//plot of biomass concentration with time
xtitle("Plot of substrate concentration", "Time, hours", "Substrate concentration, g/g");
myfile = "bioreactor-substrate.png";//file name
filename = dirbioreactor + myfile;
xs2png(1, filename);//save file as PNG

scf(2);//create new graphic window
plot(t, y3);//plot of biomass concentration with time
xtitle("Plot of product concentration", "Time, hours", "Product concentration, g/g");
myfile = "bioreactor-product.png";//file name
filename = dirbioreactor + myfile;
xs2png(2, filename);//save file as PNG