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

I am trying to run this system of ODEs for a pyrolysis reactor I am trying to model. I am unable to get an output even with using ODE15s. Please help!!

function [] = CFBR()
clear
clc
format long
zspan = .01:1:12.01;
ICs = [100,0,0,0,79,0,0,0,0,0,21];
options = optimset('display','off');
[z,F] = ode15s(@odesys,zspan,ICs,options)
C_B = F(:,1);
C_C = F(:,2);
C_G = F(:,3);
C_T = F(:,4);
C_O2 = F(:,5);
C_CO = F(:,6);
C_H2O = F(:,7);
C_CH4 = F(:,8);
C_CO2 = F(:,9);
C_H2 = F(:,10);
C_N2 = F(:,11);
plot(zspan,C_B,zspan,C_C,zspan,C_G,zspan,C_T,zspan,C_O2,zspan,C_CO,zspan,C_H2O,zspan,C_CH4,zspan,C_CO2,zspan,C_H2,zspan,C_N2)
end
function [dsdz] = odesys(z,s)
%---------- List of Thermodynamic and Kinetic Data for Reactions-----------
% A = pre-exponential factor [m/s]
% B = exponential factor [K^-1]
% H = heat of reaction [J/kg biomass]
% E = activation energy [J/kmol]
A_p11 = 1.1183e-10;
A_p12 = 8.10156e-10; 
B_p11 = 68.4868;
B_p12 = 104.5661;
deltaH_p1 = -418e3;
A_p2 = 5.9e7;
E_p2 = 123.5;
deltaH_p2 = 42e3;
A_c1 = 9.2e6;
E_c1 = 80.235;
A_c2 = 9.2e6;
E_c2 = 80.235e6;
A_c3 = 1.3e11;
E_c3 = 125.59e6;
A_c4 = 1e11;
E_c4 = 83.144e6;
A_wg = 2.78;
E_wg = 12.58e6;
A_E = .0265;
E_E = 32.9e6;
A_c5 = 5.67e9;
E_c5 = 160e6;
A_g1 = 1e7;
E_g1 = 218e6;
A_g2 = 1e5;
E_g2 = 218e6;
A_g3 = 1e5;
E_g3 = 218e6;
%----------Parameters----------
R = 8314; % Gas Constant [J/kmol-K]
dp = .05; % Biomass particle diameter [m]
av_b = 6/dp; % Specific particle surface area [m^-1]
epsilon = .5; % Overall void fraction [unit-less]
epsilon_sb = .5; % Biomass void fraction [unit-less]
rho_b = 700; % Biomass density [kg/m^3]
Ts = 550; % Temperature [K]
T = 550; % Temperature [K]
Tg = 550; % Temperature [K]
av_c = av_b;
ug = 6; % Gas velocity [m/s]
v = 5.3e-5; % Particle viscosity [Pa-s]
D = 2e-5; % Diffusivity [m^2/s]
Mc = 95; % Char molecular weight [kg/kmol]
Mt = 95; % Tar molecular weight [kg/kmol]
%----------Reaction Rate Constant Equations-----------
k_p1 = A_p11*exp(Ts/B_p11)+A_p12*exp(Ts/B_p12);
k_p2 = A_p2*exp(-E_p2/(R*T));
k_c1 = A_c1*exp(-E_c1/(R*T));
k_c2 = A_c2*exp(-E_c2/(R*T));
k_c3 = A_c3*exp(-E_c3/(R*T));
k_c4 = A_c4*exp(-E_c4/(R*T));
k_wg = A_wg*exp(-E_wg/(R*T));
K_wg = A_E*exp(-E_E/(R*T));
k_c5 = A_c5*exp(-E_c5/(R*Ts));
k_g1 = A_g1*exp(-E_g1/(R*Ts));
k_g2 = A_g2*exp(-E_g2/(R*Ts));
k_g3 = A_g3*exp(-E_g3/(R*Ts));
%----------Heat and Mass Transfer Coefficients------------
Re = (rho_b*ug*dp)/v;
Sc = v/(rho_b*D);
k_m = 2.06*ug*(Re^-.575)*(Sc^(-2/3));
%----------Defining 'dummy variables'-----------
C_B = s(1); % Biomass concentration
C_C = s(2); % Char concentration
C_G = s(3); % Gas concentration
C_T = s(4); % Tar concentration
C_O2 = s(5); % Oxygen concentration
C_CO = s(6); % CO concentration
C_H2O = s(7); % H2O concentration
C_CH4 = s(8); % CH4 concentration
C_CO2 = s(9); % CO2 concentration
C_H2 = s(10); % H2 concentration
C_N2 = s(11); % N2 concentration
%----------Reaction Rates-----------
r_p1 = k_p1*av_b*epsilon_sb*rho_b;
r_p2 = k_p2*epsilon*C_T;
r_c1 = k_c1*epsilon*Tg*(C_T^.5)*C_O2;
r_c2 = k_c2*epsilon*Tg*(C_CH4^.5)*C_O2;
r_c3 = k_c3*epsilon*C_CO*(C_O2^.25)*C_H2O;
r_c4 = k_c4*epsilon*C_H2*C_O2;
r_wg = k_wg*epsilon*(C_CO*C_H2O-((C_CO2*C_H2)/K_wg));
r_c5 = av_c*epsilon*(((k_m^-1)+(k_c5^-1))^-1)*C_O2;
r_g1 = av_c*epsilon*(((k_m^-1)+(k_g1^-1))^-1)*C_CO2;
r_g2 = av_c*epsilon*(((k_m^-1)+(k_g2^-1))^-1)*C_H2;
r_g3 = av_c*epsilon*(((k_m^-1)+(k_g3^-1))^-1)*C_H2O;
%----------Material Balances-----------
dC_B = -r_p1;
dC_C = .33*r_p1-Mc*(r_c5+r_g1+r_g2+r_g3);
dC_G = .48*r_p1+r_p2;
dC_T = .19*r_p1-r_p2-Mt*r_c1;
dC_O2 = -.871*r_c1-1.5*r_c2-r_c3-r_c4-r_c5;
dC_CO = (.156/28.01)*.48*r_p1+(Mt/28.01)*.5*r_p2+r_c1+r_c2-2*r_c3-r_wg+2*r_g1+r_g3;
dC_H2O = (.521/18.015)*.48*r_p1+.761*r_c1+2*r_c2+2*r_c4-r_wg-r_g3;
dC_CH4 = (.031/16.04)*.48*r_p1+(Mt/16.04)*.2*r_p2-r_c2+r_g2;
dC_CO2 = (.271/44.01)*.48*r_p1+(Mt/44.01)*.3*r_p2+2*r_c3+r_wg+r_c5-r_g1;
dC_H2 = (.021/2.016)*.48*r_p1-2*r_c4+r_wg-2*r_g2+r_g3;
dC_N2 = 0;
dsdz = [dC_B dC_C dC_G dC_T dC_O2 dC_CO dC_H2O dC_CH4 dC_CO2 dC_H2 dC_N2]';
end

"Warning: Failure at t=1.000000e-02. Unable to meet integration tolerances without reducing the step size below the smallest value allowed

(2.775558e-17) at time t."

[Star Strider]

Your function has encountered a singularity (+Inf or -Inf) at that point. There’s not much we can do. You need to be certain your ODE function is coded correctly.

Your code could be made more efficient, such as recoding:

r_c5 = av_c*epsilon*(((k_m^-1)+(k_c5^-1))^-1)*C_O2;

as:

r_c5 = (C_O2*av_c*epsilon*k_c5*k_m)/(k_c5 + k_m);

and similarly for the others.

[Walter Roberson]

The message can sometimes occur when there is no true singularity, in areas that the slope is steep enough. Sometimes reducing the maximum step size can help, but often that only postpones the problem: there will always be systems which have slopes that are too steep for double precision numeric methods.

Sometimes what you need to do is switch to a "stiff" solver such as ode23s.