function dy = modello(t,var)
Y1 = 1 ;
Y2 = 0.535 ;
Y3 = 1 ;
Y4 = 0.0089 ;
Y5 = 0.0152 ;
Y6 = 0.0155 ;
Y7 = 0.285 ;
Y8 = 1 ;
Y9 = 0.97 ;
Y10 = 0.99 ;
mu_max1 = 2.11 ;
mu_max2 = 0.26 ;
mu_max3 = 25.97 ;
ks1 = 4.66 ;
ks2 = 0.26 ;
ki2 = 5.7 ;
ks3 = 0.75 ;
kla1 = 0.294*24 ;
kla2 = 6.6*24 ;
kla3 = 0.261*24 ;
H_co2 = 7.48 ;
H_h2 = 0.0016 ;
H_ch4 = 0.016 ;
Vg = 9.3 ;
Vl = 52.7 ;
T = 328 ;
P = 1 ;
pw = 0.155 ;
R = 0.0821 ;
qh2_in = 1.2*24 ;
S1_in = 41.619 ;
S2_in = 0 ;
S3G_in = 0;
S5G_in = 0;
q_l = 0.448*24 ;
PM_co2 = 44 ;
PM_h2 = 2 ;
PM_ch4 = 16 ;
S4G_in = 0.074 ;
D = 1/30 ;
X1=var(1) ;
X2=var(2) ;
X3=var(3) ;
S1=var(4) ;
S2=var(5) ;
S3=var(6) ;
S4=var(7) ;
S5=var(8) ;
S3G=var(9) ;
S4G=var(10) ;
S5G=var(11) ;
mu1=mu_max1.*S1./(S1+ks1);
mu2=mu_max2.*S2./(S2+ks2+(S2^2)/ki2);
mu3=mu_max3.*S4./(S4+ks3);
S3star=(S3G./(S3G+S4G+S5G)).*P*H_co2 ;
S4star=(S4G./(S3G+S4G+S5G)).*P*H_h2 ;
S5star=(1-(S3G./(S3G+S4G+S5G))-(S4G./(S3G+S4G+S5G))).*P*H_ch4 ;
q_out=R*T/(P) *Vl*(((kla1/PM_co2)*(S3-S3star))-((kla2/PM_h2)*(S4-S4star))+(kla3/PM_ch4)*(S5-S5star));
q_in=qh2_in;
dx=zeros(11,1) ;
dx(1,1)=(mu1-D-0.05*mu_max1).*X1 ;
dx(2,1)=(mu2-D-0.027*mu_max2).*X2 ;
dx(3,1)=(mu3-D-0.05*mu_max3).*X3 ;
dx(4,1)=(D).*(S1_in - S1)-Y1.*mu1.*X1 ;
dx(5,1)=-(D).*S2 + Y2.*mu1.*X1-Y3.*mu2.*X2 ;
dx(6,1)=+Y4.*mu1.*X1+Y5.*mu2.*X2-Y6.*mu3.*X3-kla1.*(S3-S3star);
dx(7,1)=+Y7.*mu1.*X1-Y8.*mu3.*X3+kla2.*(S4star-S4);
dx(8,1)=+Y9.*mu2.*X2+Y10.*mu3.*X3-kla3.*(S5-S5star);
dx(9,1)=q_in.*S3G_in./Vg -q_out.*S3G./Vg +(Vl/Vg).*kla1.*(S3-S3star);
dx(10,1)=q_in.*S4G_in./Vg -q_out.*S4G./Vg +(Vl/Vg).*kla2.*(S4-S4star);
dx(11,1)=q_in.*S5G_in./Vg -q_out.*S5G./Vg +(Vl/Vg).*kla3.*(S5-S5star);
end
clc ; clear all; close all
range=[1:1:40];
IC=[2.06;0.68;0.01;41.619;0;0;0;0;0;0.0740;0];
[time, sol]=ode45(@Prova1,range,IC);