ode45 not solving past the initial conditions (reactor dynamics)

I first attepted to plot y(1) through y(6) but no data showed up on the figure. Then I checked the array [z, y] of the ode45 output and only the initial conditions showed up in row 1 with the rest of the array filled in with NaN. I am pretty new to ode45 but I got ode45 to solve a similar less complicated problem when I was first trying it out. Hopefully there is just some typo but I cant find it and I think this is my last option. Thanks!

zspan=[0 1];

z0=[3.1135; 0.76830; 0; 0; 0; 0.1617; 0];

[z, y]=ode45(@diffy, zspan, z0)

%y(1)=Fh2

%y(2)=Fco2

%y(3)=Fm

%y(4)=Fw

%y(5)=Fdme

%y(6)=Fco

%y(7)=H

%T=y(7)/(y(1)*cph+y(2)*cpco2+y(3)*cpm+y(4)*cpw+y(5)*cpdme+y(6)*cpco)

function dy=diffy(~,y)

%thermo kJ/mol

dh1=-56.33;

dh2=-21.255;

dh3=-40.9;

R=8.3145; %kJ/kmol/K

%reactor

eps=0.33;

rho=1900;

vdot=2.37;

Re=0.019;

Ri=0.01;

U=20;

Tv=475;

Permw=3*10^-7;

Permh=2.2*10^-7;

%Temp and cp kJ/kmol/K

cph=14.5*2;

cpco2=43.99;

cpm=87.49;

cpw=2.25*18;

cpdme=125000;

cpco=27.995*1.185;

T=y(7)/(y(1)*cph+y(2)*cpco2+y(3)*cpm+y(4)*cpw+y(5)*cpdme+y(6)*cpco);

%catalyst kinetics

k1=35.45*exp(-1.7069*10^4/R./T);

k2=8.2894*10^4*exp(-5.2940*10^4/R./T);

k3=7.3976*exp(-2.0436*10^4/R./T);

Kh=0.249*exp(3.4394*10^4/R./T);

Kco2=1.02*10^-7*exp(6.74*10^4/R./T);

Kco=7.99*10^-7*exp(5.81*10^4/R./T);

K1=exp(4213./T-5.752*log(T)-1.707*10^-3*T+2.682*10^-6*T.^2-7.232*10^-10.*T.^3+17.6);

K2=exp(4019./T+3.707*log(T)-2.783*10^-3*T+3.8*10^-7*T.^2-6.561*10^4./T.^3-26.64);

K3=10^(2167./T-0.5194*log10(T)+1.037*10^-3.*T-2.331*10^-7.*T.^2-1.2777);

%pressure conversion

Ph=y(1)./(vdot*R*T);

Pco2=y(2)./(vdot*R*T);

Pm=y(3)./(vdot*R*T);

Pw=y(4)./(vdot*R*T);

Pdme=y(5)./(vdot*R*T);

Pco=y(6)./(vdot*R*T);

%reactions

r1=k1*(Pco2*Ph*(1-Pm*Pw/K1/Pco2/Ph^3))/(1+Kco2*Pco2+Kco*Pco+sqrt(Kh*Ph))^3;

r2=k2*(Pm^2/Pw-Pdme/K2);

r3=k3*(Pw-(Pco2*Ph/K3/Pco))/(1+Kco2*Pco2+Kco*Pco+sqrt(Kh*Ph));

%perm

pw=0.05*Pw;

ph=0.05*Ph;

Jw=Permw*(Pw-pw);

Jh=Permh*(Ph-ph);

%deez

d=[0 0 0 0 0 0 0];

d(1)=-rho*(1-eps)*pi*(Re^2-Ri^2)*(3*r1-r3)-Jh*2*pi*Ri;

d(2)=-rho*(1-eps)*pi*(Re^2-Ri^2)*(r1-r3);

d(3)=rho*(1-eps)*pi*(Re^2-Ri^2)*(r1-2*r2);

d(4)=rho*(1-eps)*pi*(Re^2-Ri^2)*(r1-r3+r2)-Jw*2*pi*Ri;

d(5)=rho*(1-eps)*pi*(Re^2-Ri^2)*r2;

d(6)=-rho*(1-eps)*pi*(Re^2-Ri^2)*r3;

d(7)=(-dh1*r1-dh2*r2-dh3*r3)*rho*(1-eps)*pi*(Re^2-Ri^2)-U*(T-Tv)*2*pi*Re;

dy=d';

end

[Torsten]

Divisions by zero everywhere:

T = 0
K2 = 0
...

Answer 1

[Sulaymon Eshkabilov]

The problem is not only with T = 0 and K2 = 0, there must be some other pitfalls with the formulations. Please check the formulations and also, if necessary, you may need to adjust error tolerances with odeset options, e.g.:

OPTs = odeset('reltol', 1e-10, 'abstol', 1e-16);
[z, y]=ode45(@diffy, zspan, z0, OPTs);