Solving a system of ODEs coupled with an algebraic equation

Hi everyone,

I have tried to couple an fsolve function with a system ODEs found in:

https://www.mathworks.com/matlabcentral/answers/458083-parameter-estimation-in-ode?s_tid=prof_contriblnk

while following Star Strider ‘Igor_Moura’ function.

The system of ODEs are solved by:

t=[0.1
0.2
0.4
0.6
0.8
1
1.5
2
3
4
5
6];

c=[0.902	0.06997	0.02463	0.00218
0.8072	0.1353	0.0482	0.008192
0.6757	0.2123	0.0864	0.0289
0.5569	0.2789	0.1063	0.06233
0.4297	0.3292	0.1476	0.09756
0.3774	0.3457	0.1485	0.1255
0.2149	0.3486	0.1821	0.2526
0.141	0.3254	0.194	0.3401
0.04921	0.2445	0.1742	0.5277
0.0178	0.1728	0.1732	0.6323
0.006431	0.1091	0.1137	0.7702
0.002595	0.08301	0.08224	0.835];

theta0=[1;1;1;1;1;1];

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

fprintf(1,'\tRate Constants:\n')
for k1 = 1:length(theta)
    fprintf(1, '\t\tTheta(%d) = %8.5f\n', k1, theta(k1))
end

pHCalculation;

tv = linspace(min(t), max(t));
Cfit = kinetics5(theta, tv);
figure(1)
plot(t, c, 'p')
hold on
hlp = plot(tv, Cfit,t,ph.*0.1,'r-');
hold on
plot(t,ph.*0.1,'r-')
hold on
plot (t,pH.*0.1,'kd')
grid
xlabel('Time')
ylabel('Concentration')
legend(hlp, 'C_1(t)', 'C_2(t)', 'C_3(t)', 'C_4(t)', 'C_5(t)','pH-Mod','Location','best')


function C=kinetics(theta,t)
c0=[1;0;0;0;0];
[T,Cv]=ode45(@DifEq,t,c0);
%
    function dC=DifEq(t,c)
    dcdt=zeros(5,1);
    dcdt(1)=-theta(1).*c(1)-theta(2).*c(1);
    dcdt(2)= theta(1).*c(1)+theta(4).*c(3)-theta(3).*c(2)-theta(5).*c(2);
    dcdt(3)= theta(2).*c(1)+theta(3).*c(2)-theta(4).*c(3)+theta(6).*c(4);
    dcdt(4)= theta(5).*c(2)-theta(6).*c(4);
    dcdt(5)=theta(1).*c(2)-theta(5);
    dC=dcdt;
    end
C=Cv(:,1:4);
end


function C=kinetics5(theta,t)
c0=[1;0;0;0;0];
[T,Cv]=ode45(@DifEq,t,c0);
%
    function dC=DifEq(t,c)
    dcdt=zeros(5,1);
    dcdt(1)=-theta(1).*c(1)-theta(2).*c(1);
    dcdt(2)= theta(1).*c(1)+theta(4).*c(3)-theta(3).*c(2)-theta(5).*c(2);
    dcdt(3)= theta(2).*c(1)+theta(3).*c(2)-theta(4).*c(3)+theta(6).*c(4);
    dcdt(4)= theta(5).*c(2)-theta(6).*c(4);
    dcdt(5)=theta(1).*c(2)-theta(5);
    dC=dcdt;
    end
C=Cv;
end

And the algebraic equation is given by:

K1   = 6.24*2;
K2   = 5.68e-5*4e1;
K3   = 1.7e-3*1.1e3;
K4   = 6.55e-8;
K5   = 5.3e-8;
x0   = 7.41e-06;

t = [0.1
0.2
0.4
0.6
0.8
1
1.5
2
3
4
5
6];

A   = [0.04879
6.650241034
9.984349357
12.80945024
15.8072316
18.85920165
21.81901476
24.15458161
23.99797436
22.9179308
22.08000295
21.3980477
];

B   = [0
3.228125432
3.039568727
2.353698976
1.924485852
1.669257876
1.577120575
1.848544338
3.227063
5.416789376
7.795831728
10.27649137
];

D    = [0.04879
        6.89871101
        13.90837796
        20.90460441
        27.81415915
        34.59918462
        41.11877195
        46.61855348
        48.32509945
        47.89047761
        47.18901391
        46.2708157
        ];
    

pH      = [8.13
        7.76
        7.43
        7.45
        7.45
        7.00
        4.39
        3.61
        3.21
        2.88
        2.88
        2.64
        ];
    
for i = 1:length(t)
% BosnianKingdom;
% A = c(1);
% B = c(2);
% D = c(3);
    
      sol(i)    = fsolve(@(x)  x + 2.* A(i) - ((B(i).* K1.*x)/(x.^2 + K1.*x + K1*K2))- 2.*((B(i).*K1*K2)/(x.^2 ...
              + K1.*x + K1*K2))-((D(i).*K3.*x)/(x.^2 + K3.*x + K3*K4))- 2.*((D(i).*K3.*K4)/(x.^2 ...
              + K3.*x + K3*K4))- K5./x, x0);
end
ph = 3 - log10(sol);
plot(t,ph.*0.1,'r--')
hold on
plot (t,pH.*0.1,'kd')
xlabel('Time (sec)')
ylabel('pH')
hold off
legend('pH-Mod','pH-Exp')

The code works well when I solve ODEs seperately and manually input the ODEs-calculated values c(1), c(2) and c(3) as columns of A, B and D. I am interested in inputing c(1), c(2) and c(3) step-wise as they are being calculated from the system of ODEs. Is it possible, or should carryon the way I am doing?

Thanks for your help.