My matlab code never ends, How can I speed up matlab code?

Hello there, I am trying to create 10 ODEs related to reaction kinetic reactions.

my matlab code runs, but running time is too long that it is not over. (estimated running time is more than 10 hrs).

How can I speed up the running time, or make it more efficient?

global KCO KH2O KCO2 Ksr Kwgs K3 ksr kwgs ktotal R T
tic;

%1 Definition of parameters
    global d L dp ri dt ro void Ro A A1 A2 s m den Ps Pt0 R FCO2 ICO2 IH2 n Fsw SN2 Tf T T298 u P0 theta Da step zspan opts F0 t i j
    %Declaration of variables for outp
    % ut parameters
    global rmax theta1 theta2 theta3 theta4 theta5 DEN0


    %1.1 Membrane (For dimensionless program these parameters do not affect the
    %results
    
   
    L=3.5e-2;                 %length of the membrane, m
    ro=0.3e-2;                %outer radius of the membrane, m
    A1=2*ro*pi;               %membrane outer area per unit lenght, m
    s = A1*L;                 %reaction area, m^2 
    m=0.25;                   


%1.3 Operating conditions
    R=8.314;                  %gas constant        
    Tf=800;                   %feed temperature, C
    T=Tf+273.15;              %feed temperature, K
    T298 = T/298;
    IH2= 3 / 24500;
    SN2=100/24500;           
    

% Given Constants
% Adsorption constants, atm^-1
KCO = 7.99e-7 * exp( 58100 / R / T) ;   % Constant for reaction 1
KH2O = 4.13e-11 * exp( 104500 / R / T) ; % actual = KH2O/KH2^0.5
KCO2 = 1.02e-7 * exp( 67400 / R / T) ;   % Constant for reaction 2

% equilibrium constants, 
Ksrtop = -21664-2.8 * (T-298)- T * (-53.19 - 2.8 * log( T298 )); 
Kwgstop = 9838.1 - 6.47 * (T-298) - T * ( 10.14 - 6.47 * log( T298 )); 
K3top = -11826 - 6.75 * (T-298) - T * ( -43.06 - 6.75 * log( T298 )); 
Ksr = exp( Ksrtop / R / T) ; %Pa^-2
Kwgs = exp( Kwgstop / R / T); %1
K3 = exp( K3top / R / T) ; %Pa^-2

% Kinetic constants
ksr = 2.69e7 * exp( -109900 / R / T) ;  % reaction 1, mol/s/kg_cat/Pa
kwgs = 7.31e8 * exp( -123400 / R /T) ; % reaction 2, mol/s/kg_cat/Pa^0.5
ktotal = 4.36e2 * exp( -65200 / R / T) ; % reaction 3, mol/s/kg_cat/Pa

% Initial Condition, atm
P_CO = 1e-6;
P_CO2 = 0.25 * 101325;
P_H2 = 0.75 * 101325;
P_CH3OH = 1e-6;
P_H2O = 1e-6;
qCO = 1e-6;
qCO2 = 1e-6;
qH2 = 1e-6;
qCH3OH = 1e-6;
qH2O = 1e-6;



% Calculate DEN (common term in reaction rate equations)
DEN = (1 + KCO * P_CO + KCO2 * P_CO2) * (sqrt(P_H2) + KH2O * P_H2O);
DEN0 = (1 + KCO * 0 + KCO2 * P_CO2) * (sqrt(P_H2) + KH2O * 0);

% Calculate reaction rates
r_1 = ksr * KCO * (P_CO * P_H2^1.5 - (P_CH3OH / P_H2^0.5 / Ksr)) / DEN;
r_2 = kwgs * KCO2 * (P_CO2 * P_H2 - P_H2O * P_CO / Kwgs) / DEN;
r_3 = ktotal * KCO2 * (P_CO2 * P_H2^1.5 - P_CH3OH * P_H2O / P_H2^1.5 / K3) / DEN;

rmax = ktotal * KCO2 * (P_CO2 * P_H2^1.5 - P_CH3OH * P_H2O / P_H2^1.5 / K3) / DEN0;

% Calculate theta
    % Ф0 values mol/(s*m^2*Pa)
Pi0CH3OH = 4.97611e-05;
Pi0H2O = 1.76929e-07;
Pi0CO = 2.2748e-06;
Pi0H2 = 1.59236e-05;
Pi0CO2 = 2.1714e-06;


%Energy j/mol
ECH3OH = 15000;
EH2O = 10000;
ECO = 13100;
EH2 = 13400;
ECO2 = 12500;

%ФPhi values Ф_i =Ф_i_0*exp(-Ei/R/T)
RR = R*T;
PiCH3OH = Pi0CH3OH*exp(-ECH3OH/RR);
PiH2O = Pi0H2O*exp(-EH2O/RR);
PiCO = Pi0CO*exp(-ECO/RR);
PiH2 = Pi0H2*exp(-EH2/RR);
PiCO2 = Pi0CO2*exp(-ECO2/RR);

thetaSub = 75994*s/(IH2/60);
theta1 = PiCH3OH*thetaSub;
theta2 = PiH2O*thetaSub;
theta3 = PiCO*thetaSub;
theta4 = PiH2*thetaSub;
theta5 = PiCO2*thetaSub;

%Da number
%Da = m*rmax/F_H2_0
Da = m*rmax/IH2*60/1000;

% Calculate the derivatives of each component
dP_CH3OH = r_1 + r_3;
dP_CO2 = -(r_2 + r_3);
dP_CO = -(r_1 - r_2);
dP_H2 = -(2 * r_1 + r_2 + 3 * r_3);
dP_H2O = r_2 + r_3;



dq_CH3OH = theta1 * (P_CH3OH/75994 - qCH3OH/75994);
dq_CO2 = (theta5 / 75994) * (P_CO2 - qCO2);
dq_CO = (theta3 / 75994) * (P_CO - qCO);
dq_H2 = (theta4 / 75994) * (P_H2 - qH2);
dq_H2O = (theta2 / 75994) * (P_H2O - qH2O);

% Set up the differential equations
dPdt = [dP_CH3OH; dP_CO2; dP_CO; dP_H2; dP_H2O; dq_CH3OH; dq_CO2; dq_CO; dq_H2; dq_H2O];



% Solve the differential equations
tspan = 0:1/100:1;   % length span of integration
P0 = [1e-6, 1, 1e-6, 3, 1e-6, 1e-6, 1e-6, 1e-6, 1e-6, 1e-6];  % Initial conditions as a column vector, flow rate & atm
Pt0 = [P0(1) / sum(P0(1:5)), P0(2) / sum(P0(1:5)), P0(3) / sum(P0(1:5)), P0(4) / sum(P0(1:5)) P0(5) / sum(P0(1:5))];
%opts = odeset('RelTol',1e-5,'AbsTol',1e-8);  %Tolerance for ode45
[t, y] = ode15s(@(t, y) reactionRates(t, y), tspan, P0);%, opts);

Warning: Failure at t=7.892978e-02. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (2.220446e-16) at time t.

% Plot the results
figure;
plot(t, y(:, 1), 'r', t, y(:, 2), 'g', t, y(:, 3), 'b', t, y(:, 4), 'm', t, y(:, 5), 'c');
xlabel('Dimensionless Length');
ylabel('Flow rate');
legend('CH3OH', 'CO2', 'CO', 'H2', 'H2O');

figure;
plot(t, y(:,6), 'r', t, y(:, 7), 'g', t, y(:, 8), 'b', t, y(:, 9), 'm', t, y(:, 10), 'c');
xlabel('Dimensionless Length');
ylabel('Flow rate');
legend('CH3OH', 'CO2', 'CO', 'H2', 'H2O');

% Function to calculate the reaction rates
function dPdt = reactionRates(t, y)
global KCO KH2O KCO2 Ksr Kwgs K3 ksr kwgs ktotal R T Da rmax theta1 theta2 theta3 theta4 theta5 SN2 IH2
    %preallocate for Pdt
   
    fy = y(1) + y(2) + y(3) + y(4) + y(5);
    P_CH3OH =  101325 * y(1) / fy;
    P_CO2 = 101325 * y(2) / fy;
    P_CO = 101325 * y(3) /fy;
    P_H2 =  101325 * y(4) / fy;
    P_H2O =  101325 * y(5) / fy;

    fq = y(6) + y(7) + y(8) + y(9) + y(10) + SN2 / IH2;
    qCH3OH = 101325 * y(6) / fq;
    qCO2 = 101325 * y(7) / fq;
    qCO = 101325 * y(8) / fq;
    qH2 = 101325 * y(9) / fq;
    qH2O = 101325 * y(10) /fq;

    H2 = P_H2^0.5;
    DEN = (1 + KCO * P_CO + KCO2 * P_CO2) * (H2 + KH2O * P_H2O);
    
    r_1 = ksr * KCO * H2 * (P_CO * P_H2 - (P_CH3OH / P_H2 / Ksr)) / DEN;
    r_2 = kwgs * KCO2 * (P_CO2 * P_H2 - P_H2O * P_CO / Kwgs) / DEN;
    r_3 = ktotal * KCO2 * H2 * H2 * H2 * (P_CO2 - P_CH3OH * P_H2O / P_H2 / P_H2 / P_H2 / K3) / DEN;

    qCH3OH = ( theta1 / 75994 ) * (P_CH3OH- qCH3OH );
    qCO2 = ( theta5 / 75994 ) * (P_H2O - qH2O);
    qCO = ( theta3 / 75994 ) *(P_CO - qCO);
    qH2 = ( theta4 / 75994 ) *(P_H2 - qH2);
    qH2O = ( theta2 / 75994 ) * (P_H2O - qH2O);

    DA = Da / rmax;
    dP_CH3OH =  DA * (r_1 + r_3) - qCH3OH; 
    dP_CO2 = - DA * (r_2 + r_3) - qCO2;
    dP_CO = - DA * (r_1 - r_2) - qCO;
    dP_H2 = - DA * (2 * r_1 + r_2 + 3 * r_3) - qH2;
    dP_H2O = DA * (r_2 + r_3) - qH2O;
    dq_CH3OH = qCH3OH; 
    dq_CO2 = qCO2;
    dq_CO = qCO;
    dq_H2 = qH2;
    dq_H2O = qH2O;

    % Calculate the derivatives
    dPdt(1) = dP_CH3OH;
    dPdt(2) = dP_CO2;
    dPdt(3) = dP_CO;
    dPdt(4) = dP_H2;
    dPdt(5) = dP_H2O;
    dPdt(6) = dq_CH3OH;
    dPdt(7) = dq_CO2;
    dPdt(8) = dq_CO;
    dPdt(9) = dq_H2;
    dPdt(10) = dq_H2O;


    % Print intermediate values
    %disp('Intermediate values:')
    %disp(t)
    %fprintf('Intermediate values: t = %f, y = [%s]\n', t, num2str(y, '%f, '));
    
    
    dPdt = [dP_CH3OH; dP_CO2; dP_CO; dP_H2; dP_H2O; dq_CH3OH; dq_CO2; dq_CO; dq_H2; dq_H2O]; %dq_CH3OH; dq_CO2; dq_CO; dq_H2; dq_H2O];
    %toc;
end

[Torsten]

Use a stiff solver (in this case ode15s) (see above).

Nonetheless, check your equations since the solver has problems, e.g. PH_2 becomes negative and thus H2 complex-valued.

[Jessica]

Thanks a lot! It was helpful to use a different solver, I'll go over it.

I recorrected my equations as well. It ends when t = 0.07, but I set tspan 0 to 1. Could you check any wrong approach in my code?

[Torsten]

I can't find an obvious error.

My guess is that the problem is somewhere hidden in the chemical reaction engineering equations (either a programming error (incompatible physical units,...) or the inherent difficulty of the equations themselves).