"Busy" when running this code with ODE45

Here is the code. I have also tried running it using ODE23 instead of 45, so I don't think that's the issue (I could be wrong though, I will take any help)

tspan = [0 300];
c0 = [0.01; 0.2; 0; 55.5; 0; 0;];
[t,C] = ode45(@ca2p1_v1,tspan,c0);
ext = (c0(1) - C(:,1))/(-1);
hold on
subplot(2,1,1)
    plot(t,ext)
    legend('Extent of Epoxidation Reaction');
    xlabel('time');
    ylabel('Extent');
subplot(2,1,2)
    plot(t,C(:,1),t,C(:,2),t,C(:,3))
    legend('Concentrations vs. Time')
    xlabel('time');
    ylabel('concentration');
    legend('C8H16','H2O2','C8H16O')

function [dCdt] = ca2p1_v1(t,y)
%assign species
    C8H16 = y(1);
    H2O2 = y(2);
    C8H16O = y(3);
    H2O = y(4);
    O2 = y(5);
    CO2 = y(6);
%define rate constants
    kE = 2.96e5;
    kD = 1.5;
    kC = 1.46e11;
%define rate laws
    rE = kE * C8H16 * H2O2;
    rD = kD * H2O2;
    rC = kC * C8H16O * O2;
%define ODE for each species
    dC8H16_dt = -rE;
    dH2O2_dt = -rE - 2*rD;
    dC8H16O_dt = rE - rC;
    dH2O_dt = rE + 2*rD + 8*rC;
    dO2_dt = rD + 11.5*rC;
    dCO2_dt = 8*rC;
%derivative vector
    dCdt = [dC8H16_dt; dH2O2_dt; dC8H16O_dt; dH2O_dt; dO2_dt; dCO2_dt];
end

[Les Beckham]

What is your question? FYI: it is normal for the Matlab status bar to show "Busy" when executing code.

Answer 1

[Torsten]

Use ode15s instead of ode45 - you have a stiff system of differential equations.

Answer 2

[Walter Roberson]

If you use ode15s and set the tspan to 15.01 or less than integration succeeds and shows most concentrations on the order of 10^-3 . But set the upper bound to 15.02 and you end earlier than that with concentrations going up to 10^8

If you use ode23s then it runs to completion.

tspan = [0 300];
c0 = [0.01; 0.2; 0; 55.5; 0; 0;];
[t,C] = ode23s(@ca2p1_v1,tspan,c0);
ext = (c0(1) - C(:,1))/(-1);
hold on
tiledlayout(4,2);
nexttile();
    plot(t,ext)
    title('Extent of Epoxidation Reaction');
    xlabel('time');
    ylabel('Extent');
nexttile();
    plot(t,C(:,1));
    title('C8H16')
    xlabel('time');
    ylabel('concentration');
nexttile();
    plot(t,C(:,2));
    title('H2O2')
    xlabel('time');
    ylabel('concentration');

nexttile();
    plot(t,C(:,3));
    title('C8H160')
    xlabel('time');
    ylabel('concentration');
nexttile();
    plot(t,C(:,4));
    title('H20')
    xlabel('time');
    ylabel('concentration');
nexttile();
    plot(t,C(:,5));
    title('O2')
    xlabel('time');
    ylabel('concentration');
nexttile();
    plot(t,C(:,6));
    title('CO2')
    xlabel('time');
    ylabel('concentration');

function [dCdt] = ca2p1_v1(t,y)
%assign species
    C8H16 = y(1);
    H2O2 = y(2);
    C8H16O = y(3);
    H2O = y(4);
    O2 = y(5);
    CO2 = y(6);
%define rate constants
    kE = 2.96e5;
    kD = 1.5;
    kC = 1.46e11;
%define rate laws
    rE = kE * C8H16 * H2O2;
    rD = kD * H2O2;
    rC = kC * C8H16O * O2;
%define ODE for each species
    dC8H16_dt = -rE;
    dH2O2_dt = -rE - 2*rD;
    dC8H16O_dt = rE - rC;
    dH2O_dt = rE + 2*rD + 8*rC;
    dO2_dt = rD + 11.5*rC;
    dCO2_dt = 8*rC;
%derivative vector
    dCdt = [dC8H16_dt; dH2O2_dt; dC8H16O_dt; dH2O_dt; dO2_dt; dCO2_dt];
end