function FedBatchInit
clear, clc, format short g, format compact
tspan = [0 1.]; % Range for the independent variable 
y0 = [24.; 0; 0; 0.8; 0]; % Initial values for the dependent variables 
disp(' Variable values at the initial point ');
disp([' t    = ' num2str(tspan(1))]);
disp('           y                  dy/dt         ');
disp([y0 ODEfun(tspan(1),y0)]);
[t,y]=ode45(@ODEfun,tspan,y0);
for i=1:size(y,2)
    disp([' Solution for dependent variable y' int2str(i)]);
    disp(['              t                  y' int2str(i)]);
    disp([t y(:,i)]);
    plot(t,y(:,i));
    title([' Plot of dependent variable y' int2str(i)]);
    xlabel(' Independent variable (t)');
    ylabel([' Dependent variable y' int2str(i)]);
    pause
end
%- - - - - - - - - - - - - - - - - - - - - -
function dYfuncvecdt = ODEfun(t,Yfuncvec);
Nx = Yfuncvec(1);
Ns = Yfuncvec(2);
Np = Yfuncvec(3);
V = Yfuncvec(4);
PR = Yfuncvec(5);
mum = .3; %Kinetic constant (1/hr)
Ks = 1; %Kinetic constant (g glucose/L)
KI = 300; %Kinetic constant (g glucose/L)^2
Yxs = .4; %(g cells/g glucose)
Yxp = .15; %(g cells/g product)
S0 = 200; %Feed sustrate concentration, initialization stage (g glucose/L)
FI = .2; %Feed flow rate, initialization stage (L/hr)
S = Ns / V; %Substrate concentration (g/L)
X = Nx / V; %Cell concentration (g/L)
munet = mum * S / (Ks + S + S ^ 2 / KI); %Production rate (g product/L-hr)
P = Np / V; %Product concentration (g/L)
X0 = 0; %Feed cell concentration (g/L)
dNxdt = FI * X0 + munet * V * X; %Quantity of cells (g)
dNsdt = FI * S0 + munet * V * X / Yxs; %Quantity of glucose (g)
dNpdt = munet * V * X / Yxp; %Quntity of fermentation products (g)
dVdt = FI; %Culture volume (L)
dPRdt = 0; %Production rate (g/hr)
dYfuncvecdt = [dNxdt; dNsdt; dNpdt; dVdt; dPRdt];