function test01
L = 808e-9;
hc= 1.986e-25;
I = 116e-3;
function C=kinetics(theta,t)
c0=[0.96437;0];
[T,Cv]=ode45(@DifEq,t,c0);
function dC=DifEq(t,c)
dcdt=zeros(2,1);
dcdt(1)= -L*I*theta(1).*(1-theta(2)).*c(1);
dcdt(2)= L*I*theta(1).*theta(2).*c(1);
dC=dcdt;
end
C=Cv;
end
t=[0
66
132
198
264
330
396
462
528
594
660];
c=[0.96437 0.061
0.81894 0.1128
0.71044 0.1568
0.63977 0.1942
0.53643 0.2259
0.46207 0.2529
0.39372 0.2736
0.33552 0.2934
0.2854 0.3101
0.24209 0.3243
0.213 0.3358];
err=[0.04622 0.001
0.03544 0.001
0.02994 0.002
0.03702 0.002
0.01933 0.001
0.02637 0.001
0.03497 0.001
0.0163 0.002
0.02451 0.002
0.0064 0.001
0.01229 0.001];
theta0=[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
tv = linspace(min(t), max(t));
Cfit = kinetics(theta, tv);
figure(1)
errorbar(t, c(1:11), err(1:11), 'p')
hold on
errorbar(t, c(12:22), err(12:22), 'p')
hold on
hlp = plot(tv, Cfit);
hold off
grid
xlim([-25 700])
ylim([0 1.05])
xlabel('Time')
ylabel('Concentration')
legend(hlp, 'C_1(t)', 'C_2(t)', 'Location','N')
end