Fitting Graph Bioprocess Monod

Good evening I try to recreate a graph similar to that of an article, using kinetic equations (differential equations) with ODE45

Could someone tell me if they see my mistake? thanks

t = [0 20 40 60 80 100 120 140 160 180 200 220 240 260 280 300];
umax = 0.313;
ksp = 28.39;
ks = 47.51;
Yxs = 0.48;
Yxp = 0.50;
kis = 308.13;
kip = 299.67;
Pxmax = 83.35;
Ppmax = 107.79;
m = 0.001;
qmax = 3.69;
alfa = 1.53;
n = 1.53;
D = 0.023;
xi = 0.3660;
si = 225;
pi = 0.00;
%a(1) = x
%a(2) = s
%a(3) = p
f = @(t,a) [D*xi-a(1)*((umax*a(2)/((ks)+(a(2))+(a(2)^2/(kis)))*((1-(a(3))/(Pxmax))^alfa)))-D; 
    -(1/Yxs)*(D*xi-a(1)*((umax*a(2)/((ks)+(a(2))+(a(2)^2/(kis)))*((1-(a(3))/(Pxmax))^alfa)))-D)-(1/Yxp)*(D*(pi-a(3))+a(1)*(qmax*(a(2))/((ksp)+(a(2))+(a(2)^2/(kip)))* (1-(a(3))/(Ppmax))^n))+(m*a(1))-D*(si-a(2)); 
    D*(pi-a(3))+a(1)* (qmax*(a(2))/((ksp)+(a(2))+(a(2)^2/(kip)))* (1-(a(3))/(Ppmax))^n)];
xt0 = [0.3660,225,0.00];
[tspan,a] = ode45(f,[0 250],xt0);    
figure
plot(tspan,a(:,1),tspan,a(:,2),tspan,a(:,3))
title('Modelado Matemático: Andrews*Levenspiel')
xlabel('Tiempo en hrs') 
ylabel('Biomass X), Sustrato (S), EtOH')

Answer 1

[Alan Stevens]

I'm not convinced the equations describe the graphs! I did the following:

tspan = 0:250;
XSP0 = [0.3660,225,0.00];

[t,XSP] = ode45(@f,tspan,XSP0);  
X = XSP(:,1);
S = XSP(:,2);
P = XSP(:,3);
figure
subplot(3,1,1)
plot(t,X),grid
title('Modelado Matemático: Andrews*Levenspiel')
xlabel('Tiempo en hrs') 
ylabel('Biomass X')
subplot(3,1,2)
plot(t,S),grid
xlabel('Tiempo en hrs') 
ylabel('Sustrato (S)')
subplot(3,1,3)
plot(t,P),grid
xlabel('Tiempo en hrs') 
ylabel('EtOH')

function dXSPdt = f(~,XSP)
        
         D = 0.023;
         Xi = 0.3660;
         Si = 225;
         Pi = 0;
         X = XSP(1);
         S = XSP(2);
         P = XSP(3);
         
         dXdt = D*Xi- X*((0.313*S/(47.51+S+S^2/308.13))*(1-P/83.35)^1.53-D);
         dPdt = D*(Pi-P) + X*(3.69*S/(28.39+S+S^2/299.67))*(1-P/107.79)^1.53;
         dSdt = D*(Si-S)-dXdt/0.48-dPdt/0.50-0.001*X;
         
         dXSPdt = [dXdt; dSdt; dPdt];

end

and got the following results

Perhaps the units aren't consistent?