here is my code:
tmax = 100; % Some Arbitrary Value
tmin = 0; % Start at t = 0
n = (tmax - tmin)/h;
t = linspace(tmin, tmax, n);
Cs = zeros(n,1);
Ce = zeros(n,1);
Ces = zeros(n,1);
Cp = zeros(n,1);
%%
% k1, k2, k3 are constants in the ODEs
k1 = 2e3;
k2 = 1e-3;
k3 = 1e1;
%%
% Initial Conditions:
Cs(1) = 1; Ce(1) = 5e-5; Ces(1) = 0; CP(1) = 0;
%%
% The coupled ODEs in question:
F_Cs = @(Cs,Ce,Ces) -k1*Cs*Ce + k2*Ces;
F_Ce = @(Cs,Ce,Ces) -k1*Cs*Ce + (k2 + k3)*Ces;
F_Ces = @(Cs,Ce,Ces) k1*Cs*Ce - (k2 + k3)*Ces;
F_Cp = @(Ces) k3*Ces;
%%
for i=1:(n-1) % calculation loop
%%
% Computing slope1 for each coupled equation:
k_1 = F_Cs(Cs(i),Ce(i),Ces(i));
l_1 = F_Ce(Cs(i),Ce(i),Ces(i));
m_1 = F_Ces(Cs(i),Ce(i),Ces(i));
n_1 = F_Cp(Ces(i));
%%
% Computing slope2 for each coupled equation:
k_2 = F_Cs(Cs(i)+0.5*h*k_1, Ce(i)+0.5*h*l_1, Ces(i)+0.5*h*m_1);
l_2 = F_Ce(Cs(i)+0.5*h*k_1, Ce(i)+0.5*h*l_1, Ces(i)+0.5*h*m_1);
m_2 = F_Ces(Cs(i)+0.5*h*k_1, Ce(i)+0.5*h*l_1, Ces(i)+0.5*h*m_1);
n_2 = F_Cp(Ces(i)+0.5*h*n_1);
%%
% Computing slope3 for each coupled equation:
k_3 = F_Cs(Cs(i)+0.5*h*k_2, Ce(i)+0.5*h*l_2, Ces(i)+0.5*h*m_2);
l_3 = F_Ce(Cs(i)+0.5*h*k_2, Ce(i)+0.5*h*l_2, Ces(i)+0.5*h*m_2);
m_3 = F_Ces(Cs(i)+0.5*h*k_2, Ce(i)+0.5*h*l_2, Ces(i)+0.5*h*m_2);
n_3 = F_Cp(Ces(i)+0.5*h*n_2);
%%
% Computing slope4 for each coupled equation:
k_4 = F_Cs((Cs(i)+k_3*h),(Ce(i)+h*l_3),(Ces(i)+h*m_3));
l_4 = F_Ce((Cs(i)+k_3*h),(Ce(i)+l_3*h),(Ces(i)+h*m_3));
m_4 = F_Ces((Cs(i)+k_3*h),(Ce(i)+l_3*h),(Ces(i)+m_3*h));
n_4 = F_Cp((Ces(i)+n_3*h));
%%
% Computing the next term in the series using the 4 slopes:
Cs(i+1) = Cs(i) + (1/6)*(k_1+2*k_2+2*k_3+k_4)*h;
Ce(i+1) = Ce(i) + (1/6)*(l_1+2*l_2+2*l_3+l_4)*h;
Ces(i+1) = Ces(i) + (1/6)*(m_1+2*m_2+2*m_3+m_4)*h;
Cp(i+1) = Cp(i) + (1/6)*(n_1+2*n_2+2*n_3+n_4)*h;
end
the variables, Cs, Ce, Ces, Cp are all exploding by the 5th iteration. I don't understand why. Please help.