ODE45 code won't run

So I have this complex reaction model that I made but it just won't run. It doesn't say that there is an error necessarily, it just tries to run for a really long time. I eventually stop the the attempt and it takes me to the ode45 function script and points at a line. This is the function and code that I'm using.

clear all; close all;
clc
%Initial Concentrations
CA0 = 0.5;  %mol/L
CB0 = 0.1;  %mol/L
Qoa = 0.000025; %L/s
Qob = 0.000025; %L/s
Q1  = Qob + Qoa;

CAii1 = Qoa*CA0/Q1; CB1 = Qob*CB0/Q1; CAiiB = 0; CAiB = 0; CAi = 0; CC1 = 0; CD1 = 0;
initCon = [CAii1, CB1, CAiiB, CAi, CAiB, CC1, CD1];
V = 10*10^-3;
k2 = 0.01;
%Residence Time
tau = V/Q1;

%Volume
tauRange = [0:0.2:tau]; %s
%ODE SOLVER
[tauRange, C1] = ode45(@diffpfr, tauRange, initCon, []);

%PLOT PFR
figure (1)
plot(tauRange, C1,'Linewidth',1)


function dcdv = diffpfr(tau, C)
a = 2; b = 0.25; c = 1; d = 0.25;
k1 = 1; k1_ = 1; k2_ = 0.25;k2 = 0.25;

%Rate Equations
rAii = - a * C(1)*C(2) + b * C(3);
rB   = - a* C(1)* C(2) +  b * C(3) - c * C(4) * C(2) + d * C(7);
rAiiB = -k1 * C(4) * C(5) - k2 * C(4) * C(6) + a * C(1)*C(2) - b * C(3);
rAi  = k1 * C(4) * C(5) + k2 * C(4) * C(6) - c * C(4) * C(2) + d * C(7);
rAiB =  c * C(4) * C(2) - d * C(7) - k1_ * C(7) - k2_ * C(7);
rC = k1 * C(3) + k1_ * C(7);
rD = k2 * C(3) + k2_ * C(7);

%differential equations
dcdv = [ rAii;rB; rAiiB; rAi;  rAiB; rC; rD];
end

[Walter Roberson]

Experiment with using ode23s() as perhaps your ode is stiff.

[Najeem Afolabi]

Yeah, I tried that too but still no luck unfortunately.

Answer 1

[Cris LaPierre]

There's something with your tau. Right now it is 200. When I make it <=30 I get it to run very quickly.

If instead I use ode15s, it works as is.

%Initial Concentrations
CA0 = 0.5;  %mol/L
CB0 = 0.1;  %mol/L
Qoa = 0.000025; %L/s
Qob = 0.000025; %L/s
Q1  = Qob + Qoa;

CAii1 = Qoa*CA0/Q1; 
CB1 = Qob*CB0/Q1; 
CAiiB = 0; 
CAiB = 0; 
CAi = 0; 
CC1 = 0; 
CD1 = 0;

initCon = [CAii1, CB1, CAiiB, CAi, CAiB, CC1, CD1];

V = 10*10^-3;
k2 = 0.01;
%Residence Time
tau = 200;%V/Q1;

%Volume
tauRange = [0 tau]; %s
%ODE SOLVER  ### Switch to ode15s
[tauRange, C1] = ode15s(@diffpfr, tauRange, initCon);

%PLOT PFR
plot(tauRange, C1,'Linewidth',1)

function dcdv = diffpfr(tau, C)
a = 2; b = 0.25; c = 1; d = 0.25;
k1 = 1; k1_ = 1; k2_ = 0.25;k2 = 0.25;

%Rate Equations
rAii = - a*C(1)*C(2) + b*C(3);
rB   = - a*C(1)*C(2) +  b*C(3) - c*C(4)*C(2) + d*C(7);
rAiiB = -k1*C(4)*C(5) - k2*C(4)*C(6) + a*C(1)*C(2) - b*C(3);
rAi  = k1*C(4)*C(5) + k2*C(4)*C(6) - c*C(4)*C(2) + d*C(7);
rAiB =  c*C(4)*C(2) - d*C(7) - k1_*C(7) - k2_*C(7);
rC = k1*C(3) + k1_*C(7);
rD = k2*C(3) + k2_*C(7);

%differential equations
dcdv = [ rAii; rB; rAiiB; rAi;  rAiB; rC; rD];
end

[Cris LaPierre]

For more guidance on choosing a solver, see the Basic Solver Selection setion of the Choose and ODE Solver page.

[Najeem Afolabi]

Ahh I see. Thank you very much. I'll apply this now.