ode45 code runs indefinitely

I used this code to model ODEs f(1)-f(5) and it worked fine, I then added a sixth ODE f(6) for pressure drop and now the code runs indeffinatley. I have tried different ODE solvers including ones for stiff ODEs but get this error

Warning: Failure at t=1.056632e-20. Unable to meet

integration tolerances without reducing the step size

below the smallest value allowed (3.753912e-35) at time t.

my code is as follows:

clc
w = [0,500];
fch40 = 894065.3018;
fh2o0 = 3833657.964;
fco0 = 383170.8436;
fh20 = 1149512.531;
fco20 = 649.8427015;
P0 = 25;
f0 = [fch40,fh2o0,fco0,fh20,fco20,P0];
[w,f] = ode45(@dfdw_code2pluspd,w,f0);
plot(w,f(:,1),w,f(:,2),w,f(:,3),w,f(:,4),w,f(:,5));
xlabel('weight of catalyst');
ylabel('Flowrate');
title('F vs w'); 
legend('CH4','H2O','CO','H2','CO2');
plot(w,f(:,6));
xlabel('p or w');
ylabel('p or w');
title('pressure drop');
legend('pressure');

with dfdw_code2pluspd being:

function f = dfdw_code2pluspd(w,f)
fch40 = 894065.3018;
fh2o0 = 3833657.964;
fco0 = 383170.8436;
fh20 = 1149512.531;
fco20 = 649.8427015;
fi0 = 38082.53202;
ftot0 = fch40+fh2o0+fco0+fh20+fco20+fi0;

fch4 = f(1); 
fh2o = f(2);
fco = f(3);
fh2 = f(4);
fco2 = f(5);
fi = fi0;
ftot = fch4+fh2o+fco+fh2+fco2+fco2+fi;

P0 = 25;
P = f(6);
Pch4 = (fch4/ftot)*P;
Ph2o = (fh2o/ftot)*P;
Pco = (fco/ftot)*P;
Ph2 = (fh2/ftot)*P;
Pco2 = (fco2/ftot)*P;

A1 = 4.225e15;
A2 = 1.955e6;
A3 = 1.020e15;
E1 = 240.1;
E2 = 67.13;
E3 = 243.9;
R = 8.314;
T0 = 900;
T = T0;
Hh2o = 88.68;
Hch4 = -38.28;
Hco = -70.61;
Hh2 = -82.90;
Bh2o = 1.77e5;
Bch4 = 6.65e-4;
Bco = 8.23e-5;
Bh2 = 6.12e-9;
Dh1 = 206; %kJ/mol
Dh2 = -41.1; %kJ/mol
Dh3 = 164.9; %kJ/mol
k1 = A1*exp(-E1/(R*T));
k2 = A2*exp(-E2/(R*T));
k3 = A3*exp(-E3/(R*T));
ke1 = A1*exp(-Dh1/(R*T));
ke2 = A2*exp(-Dh2/(R*T));
ke3 = A3*exp(-Dh3/(R*T));
kch4 = Bch4*exp(-Hch4/(R*T));
kco = Bco*exp(-Hco/(R*T));
kh2 = Bh2*exp(-Hh2/(R*T));
kh2o = Bh2o*exp(-Hh2o/(R*T));

D = 0.1;
Ac = (pi*(D^2))/4;
G = ftot/Ac;
ro0 = 1.225;
roc = 1300;
mu = 4.6e-5;
Dp = 0.0015;
phi = 0.37;
beta0 = (G/(ro0*Dp))*((1-phi)/(phi^3))*(((150*(1-phi)*mu)/Dp)+(1.75*G));
alpha = (2*beta0)/(Ac*(1-phi)*roc*P0);

DEN = 1+(kch4*Pch4)+(kco*Pco)+(kh2*Ph2)+((Ph2o*kh2o)/Ph2);
r1 = (k1/(Ph2^2.5))*((((Pch4*Ph2o)-(Pco*(Ph2^3))/ke1))/(DEN^2));
r2 = (k2/(Ph2))*((((Pco*Ph2o)-(Pco2*(Ph2))/ke2))/(DEN^2));
r3 = (k3/(Ph2^3.5))*((((Pch4*(Ph2o^2))-(Pco2*(Ph2^4))/ke3))/(DEN^2));

rch4 = -r1-r3;
rh2o = -r1-r2;
rco = r1-r2;
rh2 = (3*r1)+r2+(4*r3);
rco2 = r2+r3;  

dPdw = -((alpha/2)*P0*(P0/P)*(ftot/ftot0)*(T/T0));

dfdw(1) = rch4;
dfdw(2) = rh2o;
dfdw(3) = rco;
dfdw(4) = rh2;
dfdw(5) = rco2;
dfdw(6) = dPdw;%dpdw

f = [dfdw(1), dfdw(2), dfdw(3), dfdw(4), dfdw(5), dfdw(6)]';
end

[Torsten]

Shouldn't it be

ftot = fch4+fh2o+fco+fh2+fco2+fi;

instead of

ftot = fch4+fh2o+fco+fh2+fco2+fco2+fi;

?

Of course, CO2 is important in our days, but ..

beta0 and alpha are incredibly high (in the order 1e20). You should check the formulae.

[Thomas Laverick]

Yes you are right about ftot, i will check my formulas. Thank you

Answer 1

[Sam Chak]

Hi @Thomas Laverick

The pressure drops rapidly. So, by the time sec.

If you look at the 6th equation (dPdw), you will notice that there is a division by P. When it reaches 0, singularity occurs. It's like falling into the "black hole" forever. Hope this explanation helps you understand what went wrong.