1

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.

0 Answers0