all ODE function returns NaN

Hi, I have a project where I need to model a reaction in a pack bed. The codes are below. It was running well before I added the differential equation for temperature change inside the unit (shown as dydz(7)). When I ran the code, it kept on giving NaN as solutions to the differential equation. Does this mean my function is too stiff for the ODE solver, and I need to split the code in a different function? Thanks.

zspan=0.1:0.1:10; 
n_co0= 39.757; %kmol/h
n_co20= 66.6; %kmol/h
n_h20= 200; %kmol/h
n_meoh0= 0; %kmol/h
n_h2o0= 0; %kmol/h
n0=[n_co0; n_co20; n_h20;n_meoh0;n_h2o0]; %kmol/h
nt0=sum(n0,'all'); %kmol/h
y0=n0./nt0;
y0=y0'; 
Pin= 50; %bar
Tin= 200+273.15; %K
IC= [y0 Pin Tin];
Ntube= 8000; 
Dtube= 0.05; %m
R=8.314; %kJ/kmol.K 
A= pi*(Dtube^2/4); %m^2

%%physical setting 
mw_co = 28.01; %kg/kmol
mw_co2 = 44.01; %kg/kmol
mw_h2 = 2.016; %kg.kmol
mw_meoh = 32.04; % kg/kmol
mw_h2o = 18.02; % kg.kmol 
mw= [mw_co; mw_co2; mw_h2; mw_meoh; mw_h2o];% kg.kmol 
M= n0.*mw; %kg/h
Mt= sum(M, 'all'); %kg/h
M_flowt= Mt/Ntube; %kg/h

%kinetic expression
fun=@(z,y) odefun_lovik_axial_pdrop_heattrf(z,y,M_flowt,mw_co,mw_co2,mw_h2,mw_meoh,mw_h2o,A,Dtube); 
[z,y]=ode15s(fun, zspan, IC); 

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

figure
plot(z,y(:,1),z,y(:,2),z,y(:,3),z,y(:,4),z,y(:,5))

figure
plot(z,y(:,6)); 

figure
plot(z,y(:,7)); 

function dydz=odefun_lovik_axial_pdrop_heattrf(z,y,M_flowt,mw_co,mw_co2,mw_h2,mw_meoh,mw_h2o,A,Dtube)
%%define variables
rho_cat= 1770; %kg/m^3
eps= 0.4; 
Dax= 0; 
t=0.002;%m
R=8.314;% kJ/kmol.K 
dp=0.0042;%m
rp= dp/2;%m
rr=Dtube/2; %m
ro=rr+t; %m
Tc=260+273.15; %K

%%kinetics setting
ka_ref= 0.499; 
kb_ref=6.62*10^-11; 
kc_ref= 3453.38;
kd_ref= 1.07; 
ke_ref= 1.22*10^10; 
Ea=17197;
Eb=124119;
Ec=0;
Ed=3669;
Ee=-94765;

%%equilibrium setting; 
K1eq_A= 3066; 
K2eq_A= 2071;
K1eq_B= 10.592;
K2eq_B= 2.029; 

%%kinetic expression
ka=ka_ref.*exp(Ea./(R.*y(7)));
kb=kb_ref.*exp(Eb./(R.*y(7)));
kc= kc_ref.*exp(Ec./(R.*y(7)));
kd= kd_ref.*exp(Ed./(R.*y(7)));
ke= ke_ref.*exp(Ee./(R.*y(7)));
K1eq= 10^((K1eq_A./y(7)) - K1eq_B); 
K2eq= 10^((K2eq_A./y(7)) - K2eq_B); 

%%to solve
%y(1)=yco
%y(2)=yco2
%y(3)=yh2
%y(4)=ymeoh
%y(5)=yh2o
%y(6)=P
%y(7)=Tr
1==y(1)+y(2)+y(3)+y(4)+y(5);
%%partial pressures and gas properties
p_co= y(1).*y(6); %bar
p_co2= y(2).*y(6); %bar
p_h2= y(3).*y(6); %bar
p_meoh= y(4).*y(6); %bar
p_h2o= y(5).*y(6); %bar
mu=67.2.*10^-7 + 0.21875.*10^-7.*y(7); 
c=y(6).*100./(R.*y(7)); %kmol/m^3
rho_g=c.*(y(1).*mw_co + y(2).*mw_co2 + y(3).*mw_h2+y(4).*mw_meoh+y(5).*mw_h2o); %kg/m^3
ug=M_flowt./(rho_g.*eps*A.*3600); %m/s
G=M_flowt./(A.*3600); %kg/m^2.s
Re=G.*dp./mu; %Reynold's number 

%%thermal properties 
ks= 4.18e-3; %...thermal conductivity
kw= 20e-3; %...thermal conductivity
Rec= 1923; %cooling water Re
Prc= 0.843; %cooling water Pr
kg= 0.01234e-3 + (1.84375e-7.*y(7)); %gas thermal conductivity
cp_co= 30.87 + (-1.285e-2.*y(7)) + (2.789e-5.*y(7).^2) + (-1.272e-8.*y(7).^3); %j/mol.K
cp_co2= 19.8 + (7.344e-2.*y(7)) + (-5.602e-5.*y(7).^2) + (1.715e-8.*y(7).^3); %j/mol.K
cp_h2= 27.14 + (9.274e-3.*y(7)) + (-1.381e-5.*y(7).^2) + (7.645e-9.*y(7).^3); %j/mol.K
cp_meoh= 21.15 + (7.092e-2.*y(7)) + (2.587e-5.*y(7).^2) + (-2.852e-8.*y(7).^3); %j/mol.K
cp_h2o= 32.24 + (1.924e-3.*y(7)) + (1.055e-5.*y(7).^2) + (3.596e-9.*y(7).^3); %j/mol.K
cpg= (y(1).*cp_co + y(2).*cp_co2 + y(3).*cp_h2 + y(4).*cp_meoh + y(5).*cp_h2o)./((y(1).*mw_co + y(2).*mw_co2 + y(3).*mw_h2 + y(4).*mw_meoh + y(5).*mw_h2o)); 
Pr= mu.*cpg./kg; %prandtl number 
gamma= ks./kg; %ratio ks:kg
kstat= kg.*((0.5 + 0.493.*((gamma-1)./(22+gamma))).*(gamma./(gamma.*eps + (1-eps))) + 1 - (0.5 + 0.493.*((gamma-1)./(22+gamma))).*(eps+(1-eps).*gamma)); 
kdyn= (kg.*Re.*Pr)./((7 + 135.8.*(rp/rr)).^2); 
krg= kdyn + kstat; %radial thermal conductivity
krs= kstat; 
hig= (1.23.*Re.*0.53.*kg)./(2.*rp); %gas inside heat transfer coefficient
his= (0.48 + 0.192.*((rr/rp)-1).^2.*krs)./rp; % inside surface heat transfer coefficient
hoc = (0.020.*Rec.^0.8.*Prc.^(1/3).*(ro/rr).^0.53.*kw)./(2.*rr); %outside cooling water heat transfer coefficient 
Uwg= ((1./hig)+((rr/ro)./hoc)+(rr.*log(rr./ro)./(2.*kw))).^-1; %overall heat transfer coefficient of gas
Uws= ((1./his)+((rr/ro)./hoc)+(rr.*log(rr./ro)./(2.*kw))).^-1; %overall heat transfer coefficient of surface
Ueff= ((1./Uwg).*(3.*rr./(8.*krg))).^-1 + ((1./Uws).*(3.*rr./(8.*krs))).^-1; %effective overall heat transfer coefficient of surface
Hr1= 5798 + 35.*(y(7)-498.15); 
Hr2= -39892 + 8.*(y(7)-498.15); 

%kinetic expression
r1=(kd.*p_co2.*p_h2.*(1-((1./K1eq).*(p_h2o.*p_meoh./(p_h2.^3.*p_co2)))))./(1+kc.*(p_h2o./p_h2)+ka.*p_h2.^0.5+kb.*p_h2o).^3;
r2=(ke.*p_co2.*((1-K2eq.*(p_h2o.*p_co./(p_co2.*p_h2)))))./(1+kc.*(p_h2o./p_h2)+ka.*(p_h2.^0.5)+kb.*p_h2o); 
 
%%differential eq
dydz(1)= ((1-eps).*rho_cat.*(r2 - -2.*r1.*y(1)))./((c.*ug.*eps).*(1-Dax.*c)); 
dydz(2)= ((1-eps).*rho_cat.*(-1.*r2 + -1.*r1 - -2.*r1.*y(2)))./((c.*ug.*eps).*(1-Dax.*c)); 
dydz(3)= ((1-eps).*rho_cat.*(-1.*r2 + -3.*r1 - -2.*r1.*y(3)))./((c.*ug.*eps).*(1-Dax.*c)); 
dydz(4)= ((1-eps).*rho_cat.*(r1 - -2.*r1.*y(4)))./((c.*ug.*eps).*(1-Dax.*c)); 
dydz(5)= ((1-eps).*rho_cat.*(r2+r1 - -2.*r1.*y(5)))./((c.*ug.*eps).*(1-Dax.*c)); 
dydz(6)= -10^-5.*(1.75+(150.*(1-eps)./Re)).*((1-eps)./eps.^3).*(ug.^2.*rho_g./dp); 
dydz(7)= (1./(G.*cpg)).*((2.*Ueff./rr).*(Tc-y(7))+((1-eps).*rho_cat.*(Hr1+Hr2)));  
dydz=dydz'; 

end

[Torsten]

Check whether the parameters used for defining dT/dt are physically senseful.

E.g. Ueff is negative at time t=0, the expression (1-eps).*rho_cat.*(Hr1+Hr2) cannot be correct because Hr1 and Hr2 must be multiplied by reaction rates etc. etc.

[Sam Chak]

First thing to check is to look for any division-by-zero terms that involve y(7), the newly added state variable. Multiples are found in code. The prime suspect is found.

Second, reduce the time step in the solver so that you see how y(7) evolves over time before the singularity occurs. Sulaymon has shown that y(7), the temperature drops so fast to nearly 0 K in a fraction of 0.1 sec.

Third, the rapid temperature drop must be caused by dT/dt or dydz(7). Some coefficients must have very large values. Are they logical for any physical processes on Earth?

[Jonathan Hanandhito]

Hi Sam and Torsten, I've managed to fix it now based on what you all have said. Some expressions were left out, and the units for some parameters and equations were wrong. Thank you so much for pointing these out!

Answer 1

[Sulaymon Eshkabilov]

There are two ways to address this issue. (1) Let ode15s to select appropriate step size; (2) make the step size much smaller than it was specified. Also, note that the solutions are within the range of [0.1 ... 0.10002]. Here is the option 1:

zspan=[0.1, 0.10002]; % Alt. way: zspan=[0.1:1e-13:.1002]; 
n_co0= 39.757; %kmol/h
n_co20= 66.6; %kmol/h
n_h20= 200; %kmol/h
n_meoh0= 0; %kmol/h
n_h2o0= 0; %kmol/h
n0=[n_co0; n_co20; n_h20;n_meoh0;n_h2o0]; %kmol/h
nt0=sum(n0,'all'); %kmol/h
y0=n0./nt0;
y0=y0'; 
Pin= 50; %bar
Tin= 200+273.15; %K
IC= [y0 Pin Tin];
Ntube= 8000; 
Dtube= 0.05; %m
R=8.314; %kJ/kmol.K 
A= pi*(Dtube^2/4); %m^2

%%physical setting 
mw_co = 28.01; %kg/kmol
mw_co2 = 44.01; %kg/kmol
mw_h2 = 2.016; %kg.kmol
mw_meoh = 32.04; % kg/kmol
mw_h2o = 18.02; % kg.kmol 
mw= [mw_co; mw_co2; mw_h2; mw_meoh; mw_h2o];% kg.kmol 
M= n0.*mw; %kg/h
Mt= sum(M, 'all'); %kg/h
M_flowt= Mt/Ntube; %kg/h

%kinetic expression
fun=@(z,y) odefun_lovik_axial_pdrop_heattrf(z,y,M_flowt,mw_co,mw_co2,mw_h2,mw_meoh,mw_h2o,A,Dtube); 
[z,y]=ode15s(fun, zspan, IC); 

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

figure
plot(z,y(:,1),z,y(:,2),z,y(:,3),z,y(:,4),z,y(:,5))

figure
plot(z,y(:,6)); 

figure
plot(z,y(:,7)); 

function dydz=odefun_lovik_axial_pdrop_heattrf(z,y,M_flowt,mw_co,mw_co2,mw_h2,mw_meoh,mw_h2o,A,Dtube)
%%define variables
rho_cat= 1770; %kg/m^3
eps= 0.4; 
Dax= 0; 
t=0.002;%m
R=8.314;% kJ/kmol.K 
dp=0.0042;%m
rp= dp/2;%m
rr=Dtube/2; %m
ro=rr+t; %m
Tc=260+273.15; %K

%%kinetics setting
ka_ref= 0.499; 
kb_ref=6.62*10^-11; 
kc_ref= 3453.38;
kd_ref= 1.07; 
ke_ref= 1.22*10^10; 
Ea=17197;
Eb=124119;
Ec=0;
Ed=3669;
Ee=-94765;

%%equilibrium setting; 
K1eq_A= 3066; 
K2eq_A= 2071;
K1eq_B= 10.592;
K2eq_B= 2.029; 

%%kinetic expression
ka=ka_ref.*exp(Ea./(R.*y(7)));
kb=kb_ref.*exp(Eb./(R.*y(7)));
kc= kc_ref.*exp(Ec./(R.*y(7)));
kd= kd_ref.*exp(Ed./(R.*y(7)));
ke= ke_ref.*exp(Ee./(R.*y(7)));
K1eq= 10^((K1eq_A./y(7)) - K1eq_B); 
K2eq= 10^((K2eq_A./y(7)) - K2eq_B); 

%%to solve
%y(1)=yco
%y(2)=yco2
%y(3)=yh2
%y(4)=ymeoh
%y(5)=yh2o
%y(6)=P
%y(7)=Tr
1==y(1)+y(2)+y(3)+y(4)+y(5);
%%partial pressures and gas properties
p_co= y(1).*y(6); %bar
p_co2= y(2).*y(6); %bar
p_h2= y(3).*y(6); %bar
p_meoh= y(4).*y(6); %bar
p_h2o= y(5).*y(6); %bar
mu=67.2.*10^-7 + 0.21875.*10^-7.*y(7); 
c=y(6).*100./(R.*y(7)); %kmol/m^3
rho_g=c.*(y(1).*mw_co + y(2).*mw_co2 + y(3).*mw_h2+y(4).*mw_meoh+y(5).*mw_h2o); %kg/m^3
ug=M_flowt./(rho_g.*eps*A.*3600); %m/s
G=M_flowt./(A.*3600); %kg/m^2.s
Re=G.*dp./mu; %Reynold's number 

%%thermal properties 
ks= 4.18e-3; %...thermal conductivity
kw= 20e-3; %...thermal conductivity
Rec= 1923; %cooling water Re
Prc= 0.843; %cooling water Pr
kg= 0.01234e-3 + (1.84375e-7.*y(7)); %gas thermal conductivity
cp_co= 30.87 + (-1.285e-2.*y(7)) + (2.789e-5.*y(7).^2) + (-1.272e-8.*y(7).^3); %j/mol.K
cp_co2= 19.8 + (7.344e-2.*y(7)) + (-5.602e-5.*y(7).^2) + (1.715e-8.*y(7).^3); %j/mol.K
cp_h2= 27.14 + (9.274e-3.*y(7)) + (-1.381e-5.*y(7).^2) + (7.645e-9.*y(7).^3); %j/mol.K
cp_meoh= 21.15 + (7.092e-2.*y(7)) + (2.587e-5.*y(7).^2) + (-2.852e-8.*y(7).^3); %j/mol.K
cp_h2o= 32.24 + (1.924e-3.*y(7)) + (1.055e-5.*y(7).^2) + (3.596e-9.*y(7).^3); %j/mol.K
cpg= (y(1).*cp_co + y(2).*cp_co2 + y(3).*cp_h2 + y(4).*cp_meoh + y(5).*cp_h2o)./((y(1).*mw_co + y(2).*mw_co2 + y(3).*mw_h2 + y(4).*mw_meoh + y(5).*mw_h2o)); 
Pr= mu.*cpg./kg; %prandtl number 
gamma= ks./kg; %ratio ks:kg
kstat= kg.*((0.5 + 0.493.*((gamma-1)./(22+gamma))).*(gamma./(gamma.*eps + (1-eps))) + 1 - (0.5 + 0.493.*((gamma-1)./(22+gamma))).*(eps+(1-eps).*gamma)); 
kdyn= (kg.*Re.*Pr)./((7 + 135.8.*(rp/rr)).^2); 
krg= kdyn + kstat; %radial thermal conductivity
krs= kstat; 
hig= (1.23.*Re.*0.53.*kg)./(2.*rp); %gas inside heat transfer coefficient
his= (0.48 + 0.192.*((rr/rp)-1).^2.*krs)./rp; % inside surface heat transfer coefficient
hoc = (0.020.*Rec.^0.8.*Prc.^(1/3).*(ro/rr).^0.53.*kw)./(2.*rr); %outside cooling water heat transfer coefficient 
Uwg= ((1./hig)+((rr/ro)./hoc)+(rr.*log(rr./ro)./(2.*kw))).^-1; %overall heat transfer coefficient of gas
Uws= ((1./his)+((rr/ro)./hoc)+(rr.*log(rr./ro)./(2.*kw))).^-1; %overall heat transfer coefficient of surface
Ueff= ((1./Uwg).*(3.*rr./(8.*krg))).^-1 + ((1./Uws).*(3.*rr./(8.*krs))).^-1; %effective overall heat transfer coefficient of surface
Hr1= 5798 + 35.*(y(7)-498.15); 
Hr2= -39892 + 8.*(y(7)-498.15); 

%kinetic expression
r1=(kd.*p_co2.*p_h2.*(1-((1./K1eq).*(p_h2o.*p_meoh./(p_h2.^3.*p_co2)))))./(1+kc.*(p_h2o./p_h2)+ka.*p_h2.^0.5+kb.*p_h2o).^3;
r2=(ke.*p_co2.*((1-K2eq.*(p_h2o.*p_co./(p_co2.*p_h2)))))./(1+kc.*(p_h2o./p_h2)+ka.*(p_h2.^0.5)+kb.*p_h2o); 
 
%%differential eq
dydz(1)= ((1-eps).*rho_cat.*(r2 - -2.*r1.*y(1)))./((c.*ug.*eps).*(1-Dax.*c)); 
dydz(2)= ((1-eps).*rho_cat.*(-1.*r2 + -1.*r1 - -2.*r1.*y(2)))./((c.*ug.*eps).*(1-Dax.*c)); 
dydz(3)= ((1-eps).*rho_cat.*(-1.*r2 + -3.*r1 - -2.*r1.*y(3)))./((c.*ug.*eps).*(1-Dax.*c)); 
dydz(4)= ((1-eps).*rho_cat.*(r1 - -2.*r1.*y(4)))./((c.*ug.*eps).*(1-Dax.*c)); 
dydz(5)= ((1-eps).*rho_cat.*(r2+r1 - -2.*r1.*y(5)))./((c.*ug.*eps).*(1-Dax.*c)); 
dydz(6)= -10^-5.*(1.75+(150.*(1-eps)./Re)).*((1-eps)./eps.^3).*(ug.^2.*rho_g./dp); 
dydz(7)= (1./(G.*cpg)).*((2.*Ueff./rr).*(Tc-y(7))+((1-eps).*rho_cat.*(Hr1+Hr2)));  
dydz=dydz'; 
end

Answer 2

[Sulaymon Eshkabilov]

Note that the found numerical solutions of y(:,1), .. y(:,5) are in different scales and therefore, it is appropriate to polt them in separate of using yyaxis option, e.g.:

zspan=[0.1,.1002]; 
%zspan=[0.1:1e-13:.1002]; 
n_co0= 39.757; %kmol/h
n_co20= 66.6; %kmol/h
n_h20= 200; %kmol/h
n_meoh0= 0; %kmol/h
n_h2o0= 0; %kmol/h
n0=[n_co0; n_co20; n_h20;n_meoh0;n_h2o0]; %kmol/h
nt0=sum(n0,'all'); %kmol/h
y0=n0./nt0;
y0=y0'; 
Pin= 50; %bar
Tin= 200+273.15; %K
IC= [y0 Pin Tin];
Ntube= 8000; 
Dtube= 0.05; %m
R=8.314; %kJ/kmol.K 
A= pi*(Dtube^2/4); %m^2

%%physical setting 
mw_co = 28.01; %kg/kmol
mw_co2 = 44.01; %kg/kmol
mw_h2 = 2.016; %kg.kmol
mw_meoh = 32.04; % kg/kmol
mw_h2o = 18.02; % kg.kmol 
mw= [mw_co; mw_co2; mw_h2; mw_meoh; mw_h2o];% kg.kmol 
M= n0.*mw; %kg/h
Mt= sum(M, 'all'); %kg/h
M_flowt= Mt/Ntube; %kg/h

%kinetic expression
fun=@(z,y) odefun_lovik_axial_pdrop_heattrf(z,y,M_flowt,mw_co,mw_co2,mw_h2,mw_meoh,mw_h2o,A,Dtube); 
[z,y]=ode15s(fun, zspan, IC); 

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

figure
subplot(311)
yyaxis left
plot(z,y(:,1))
ylabel('y(:,1)')
yyaxis right
plot(z,y(:,2))
ylabel('y(:,2)')
subplot(312)
yyaxis left
plot(z,y(:,3))
ylabel('y(:,3)')
yyaxis right
plot(z,y(:,4))
ylabel('y(:,4)')
subplot(313)
plot(z,y(:,5))
ylabel('y(:,5)')

figure
plot(z,y(:,6)); 

figure
plot(z,y(:,7)); 

function dydz=odefun_lovik_axial_pdrop_heattrf(z,y,M_flowt,mw_co,mw_co2,mw_h2,mw_meoh,mw_h2o,A,Dtube)
%%define variables
rho_cat= 1770; %kg/m^3
eps= 0.4; 
Dax= 0; 
t=0.002;%m
R=8.314;% kJ/kmol.K 
dp=0.0042;%m
rp= dp/2;%m
rr=Dtube/2; %m
ro=rr+t; %m
Tc=260+273.15; %K

%%kinetics setting
ka_ref= 0.499; 
kb_ref=6.62*10^-11; 
kc_ref= 3453.38;
kd_ref= 1.07; 
ke_ref= 1.22*10^10; 
Ea=17197;
Eb=124119;
Ec=0;
Ed=3669;
Ee=-94765;

%%equilibrium setting; 
K1eq_A= 3066; 
K2eq_A= 2071;
K1eq_B= 10.592;
K2eq_B= 2.029; 

%%kinetic expression
ka=ka_ref.*exp(Ea./(R.*y(7)));
kb=kb_ref.*exp(Eb./(R.*y(7)));
kc= kc_ref.*exp(Ec./(R.*y(7)));
kd= kd_ref.*exp(Ed./(R.*y(7)));
ke= ke_ref.*exp(Ee./(R.*y(7)));
K1eq= 10^((K1eq_A./y(7)) - K1eq_B); 
K2eq= 10^((K2eq_A./y(7)) - K2eq_B); 

%%to solve
%y(1)=yco
%y(2)=yco2
%y(3)=yh2
%y(4)=ymeoh
%y(5)=yh2o
%y(6)=P
%y(7)=Tr
%1==y(1)+y(2)+y(3)+y(4)+y(5);
%%partial pressures and gas properties
p_co= y(1).*y(6); %bar
p_co2= y(2).*y(6); %bar
p_h2= y(3).*y(6); %bar
p_meoh= y(4).*y(6); %bar
p_h2o= y(5).*y(6); %bar
mu=67.2.*10^-7 + 0.21875.*10^-7.*y(7); 
c=y(6).*100./(R.*y(7)); %kmol/m^3
rho_g=c.*(y(1).*mw_co + y(2).*mw_co2 + y(3).*mw_h2+y(4).*mw_meoh+y(5).*mw_h2o); %kg/m^3
ug=M_flowt./(rho_g.*eps*A.*3600); %m/s
G=M_flowt./(A.*3600); %kg/m^2.s
Re=G.*dp./mu; %Reynold's number 

%%thermal properties 
ks= 4.18e-3; %...thermal conductivity
kw= 20e-3; %...thermal conductivity
Rec= 1923; %cooling water Re
Prc= 0.843; %cooling water Pr
kg= 0.01234e-3 + (1.84375e-7.*y(7)); %gas thermal conductivity
cp_co= 30.87 + (-1.285e-2.*y(7)) + (2.789e-5.*y(7).^2) + (-1.272e-8.*y(7).^3); %j/mol.K
cp_co2= 19.8 + (7.344e-2.*y(7)) + (-5.602e-5.*y(7).^2) + (1.715e-8.*y(7).^3); %j/mol.K
cp_h2= 27.14 + (9.274e-3.*y(7)) + (-1.381e-5.*y(7).^2) + (7.645e-9.*y(7).^3); %j/mol.K
cp_meoh= 21.15 + (7.092e-2.*y(7)) + (2.587e-5.*y(7).^2) + (-2.852e-8.*y(7).^3); %j/mol.K
cp_h2o= 32.24 + (1.924e-3.*y(7)) + (1.055e-5.*y(7).^2) + (3.596e-9.*y(7).^3); %j/mol.K
cpg= (y(1).*cp_co + y(2).*cp_co2 + y(3).*cp_h2 + y(4).*cp_meoh + y(5).*cp_h2o)./((y(1).*mw_co + y(2).*mw_co2 + y(3).*mw_h2 + y(4).*mw_meoh + y(5).*mw_h2o)); 
Pr= mu.*cpg./kg; %prandtl number 
gamma= ks./kg; %ratio ks:kg
kstat= kg.*((0.5 + 0.493.*((gamma-1)./(22+gamma))).*(gamma./(gamma.*eps + (1-eps))) + 1 - (0.5 + 0.493.*((gamma-1)./(22+gamma))).*(eps+(1-eps).*gamma)); 
kdyn= (kg.*Re.*Pr)./((7 + 135.8.*(rp/rr)).^2); 
krg= kdyn + kstat; %radial thermal conductivity
krs= kstat; 
hig= (1.23.*Re.*0.53.*kg)./(2.*rp); %gas inside heat transfer coefficient
his= (0.48 + 0.192.*((rr/rp)-1).^2.*krs)./rp; % inside surface heat transfer coefficient
hoc = (0.020.*Rec.^0.8.*Prc.^(1/3).*(ro/rr).^0.53.*kw)./(2.*rr); %outside cooling water heat transfer coefficient 
Uwg= ((1./hig)+((rr/ro)./hoc)+(rr.*log(rr./ro)./(2.*kw))).^-1; %overall heat transfer coefficient of gas
Uws= ((1./his)+((rr/ro)./hoc)+(rr.*log(rr./ro)./(2.*kw))).^-1; %overall heat transfer coefficient of surface
Ueff= ((1./Uwg).*(3.*rr./(8.*krg))).^-1 + ((1./Uws).*(3.*rr./(8.*krs))).^-1; %effective overall heat transfer coefficient of surface
Hr1= 5798 + 35.*(y(7)-498.15); 
Hr2= -39892 + 8.*(y(7)-498.15); 

%kinetic expression
r1=(kd.*p_co2.*p_h2.*(1-((1./K1eq).*(p_h2o.*p_meoh./(p_h2.^3.*p_co2)))))./(1+kc.*(p_h2o./p_h2)+ka.*p_h2.^0.5+kb.*p_h2o).^3;
r2=(ke.*p_co2.*((1-K2eq.*(p_h2o.*p_co./(p_co2.*p_h2)))))./(1+kc.*(p_h2o./p_h2)+ka.*(p_h2.^0.5)+kb.*p_h2o); 
 
%%differential eq
dydz(1)= ((1-eps).*rho_cat.*(r2 - -2.*r1.*y(1)))./((c.*ug.*eps).*(1-Dax.*c)); 
dydz(2)= ((1-eps).*rho_cat.*(-1.*r2 + -1.*r1 - -2.*r1.*y(2)))./((c.*ug.*eps).*(1-Dax.*c)); 
dydz(3)= ((1-eps).*rho_cat.*(-1.*r2 + -3.*r1 - -2.*r1.*y(3)))./((c.*ug.*eps).*(1-Dax.*c)); 
dydz(4)= ((1-eps).*rho_cat.*(r1 - -2.*r1.*y(4)))./((c.*ug.*eps).*(1-Dax.*c)); 
dydz(5)= ((1-eps).*rho_cat.*(r2+r1 - -2.*r1.*y(5)))./((c.*ug.*eps).*(1-Dax.*c)); 
dydz(6)= -10^-5.*(1.75+(150.*(1-eps)./Re)).*((1-eps)./eps.^3).*(ug.^2.*rho_g./dp); 
dydz(7)= (1./(G.*cpg)).*((2.*Ueff./rr).*(Tc-y(7))+((1-eps).*rho_cat.*(Hr1+Hr2)));  
dydz=dydz'; 

end