how do i solve six ODE and optimization

Hello team math works,

How can i solve ode system contain 6 ode then optimize to 13 parameters. i solved this system by using Igor_moura but the results was wrong.

the model is

X1,X2, S1,S2, C,Z are variables

k1, k2, k3, k4 k5 ,k6 ,kI2, Ks1, Ks2, kla , alpha, meu1max, meu2max are parameters which optimized

D, kh, kb ,s1in, s2in, cin, zin , pt are constants

the code is below

function modelacion

function C=kinetics(theta,t)

c0=[315;2.95;1;0.92;220;100];

[T,Cv]=ode45(@DifEq,t,c0);

function dC=DifEq(t,c)

zin = 579.56;

sIin = 73.14;

sIIin = 379.8;

cin = 200.56;

D = 0.8;

kh = 23;

kb = 0.00065;

pt = 1;

dcdt=zeros(6,1);

dcdt(1)= D*(zin-c(1));

dcdt(2)= ((theta(2)*(c(3)/(c(3)+theta(4))))-(theta(1)*D))*c(2);

dcdt(3)= D*(sIin-c(3))-theta(6)*(theta(2)*(c(3)/(c(3)+theta(4))))*c(2);

dcdt(4)= ((theta(3)*(c(5)/(c(5)+theta(5)+((c(5)).^2/theta(12)))))-(theta(1)*D))*c(4);

dcdt(5) = D*(sIIin-c(5))-theta(7)*(theta(2)*(c(3)/(c(3)+theta(4))))*c(2)-theta(8)*(theta(3)*(c(5)/(c(5)+theta(5)+((c(5)).^2/theta(12)))))*c(4);

dcdt(6) = D*(cin-c(6))-theta(13)*(c(6)+c(5)-c(1)-(kh*((c(6)+c(5)-c(1)+(kh*pt)+(theta(11)/theta(13))*(theta(3)*(c(5)/(c(5)+theta(5)+(((c(5)).^2)/theta(12)))))*c(4))-(sqrt(((c(6)+c(5)-c(1)+(kh*pt)+(theta(11)/theta(13))*(theta(3)*(c(5)/(c(5)+theta(5)+(((c(5)).^2)/theta(12)))))*c(4)).^2)-(4*kh*pt*(c(6)+c(5)-c(1))))))/(2*kh)))+theta(7)*(theta(2)*(c(3)/(c(3)+theta(4))))*c(2)+theta(10)*(theta(3)*(c(5)/(c(5)+theta(5)+(((c(5)).^2)/theta(12)))))*c(4);

dC=dcdt;

end

C=Cv;

end

t=[0;1;2;3;4;5;6;7;8;9;10;11;12;13;14;15;16;17;18;19;20;21;22;23;24;25;26;27;28;29;30;31;32;33;34;35;36];

c = [315 2.95 1 0.92 220 100

314 2.88 0.5 0.89 214 103

317 2.86 0.8 0.87 214 109

320 2.8 1 0.86 214 111

320 2.78 1.1 0.85 214 112

320 2.8 1.1 0.86 214 113

317 2.78 0.8 0.86 205 114

315 2.8 0.6 0.84 198 120

317 2.82 1 0.83 200 124

319 2.8 1.1 0.82 200 124

320 2.8 1 0.83 200 124

319 2.78 1.1 0.83 205 124

319 2.76 1.1 0.84 205 122

315 2.75 0.9 0.85 198 122

313 2.79 0.5 0.84 190 129

317 2.8 0.9 0.83 191 130

317 2.78 1.1 0.81 198 130

315 2.76 0.7 0.82 190 131

316 2.76 0.8 0.8 192 131

317 2.76 1 0.8 195 131

315 2.77 0.8 0.82 190 131

315 2.76 0.6 0.8 183 136

313 2.78 1 0.78 182 137

317 2.75 1.1 0.76 195 135

318 2.69 1.2 0.76 190 132

319 2.68 1.6 0.79 200 130

319 2.68 1.1 0.81 205 125

316 2.7 1 0.83 196 123

317 2.78 0.6 0.82 195 126

318 2.8 0.7 0.81 195 130

318 2.82 1 0.81 195 130

315 2.82 0.8 0.83 190 130

317 2.8 1.1 0.8 192 130

316 2.79 1.1 0.8 200 129

315 2.8 0.8 0.83 192 130

312 2.78 0.7 0.8 185 131

320 2.7 1.2 0.77 196 131];

theta0=[1.2;0.8;0.5;9;33.7;26;226.6;637.6;34.9;24.8;155;250;80.4];

[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c);

fprintf(1,'\tConstants:\n')

for k1 = 1:length(theta)

fprintf(1, '\t\tTheta(%d) = %8.5f\n', k1, theta(k1))

end

tv = linspace(min(t), max(t));

Cfit = kinetics(theta, tv);

figure(1)

plot(t,c (:,1),'p')

hold on

hlp(1)=plot(tv,Cfit(:,1));

hold off

ylabel('Concentration (g/L)')

figure(2)

plot(t,c (:,2),'p')

hold on

hlp(2)=plot(tv,Cfit(:,2));

hold off

ylabel('Concentration (g/L)')

figure(3)

plot(t,c (:,3),'p')

hold on

hlp(3)=plot(tv,Cfit(:,3));

hold off

ylabel('Concentration (g/L)')

figure(4)

plot(t,c (:,4),'p')

hold on

hlp(4)=plot(tv,Cfit(:,4));

hold off

ylabel('Concentration (g/L)')

figure(5)

plot(t,c (:,5),'p')

hold on

hlp(5)=plot(tv,Cfit(:,5));

hold off

ylabel('Concentration (g/L)')

figure(6)

plot(t,c (:,6),'p')

hold on

hlp(6)=plot(tv,Cfit(:,6));

hold off

ylabel('Concentration (g/L)')

grid

xlabel('Time')

ylabel('Concentration')

legend(hlp, 'C_1(t)', 'C_2(t)', 'C_3(t)','C_4(t)','C_5(t)','C_6(t)', 'Location','N')

end