Why ode45 can not solve while ode15s can solve?

My code can solve with ode15s but I would like to try on ode45, So can I use RelTol with value 1e-8 or can you give me other value of RelTol?

P0                                    =  [0.5]                                                                                                                        ;%bar         (Intial pressure)
Ta0                                   =  [873]                                                                                                                  ;%K           (Intial air heat temperature)              
repeat                                =  [1]                                                                                                                          ;%Repeat condition Ta constant and Ta non constant
sheet                                 =  0                                                                                                                              ;%Sheet in excel 
condition                             =  0                                                                                                                              ;%Number of condition
for a = 1:length(repeat)  
    for i = 1:length(P0)
        for j = 1:length(Ta0)
            T0                            = Ta0(j)                                                                                                                      ;%K           (Intial temperature before entering bed)
            pj0                           = yj0.*P0(i)                                                                                                                  ;%bar         (Partial pressure)
            Fj0                           = (pj0.*v0_feed)./(R_feed.*T0)                                                                                                ;%mol/s       (Molar flowrate)
            Lspan                         = [0:0.01:0.03]                                                                                                              ;%m           (Bed length range)
            dep_var0                      = [Fj0' T0 Ta0(j)]                                                                                                            ;%[mol/s,K,K] (Intial dependent variable)
            condition                     = condition + 1                                                                                                               ;
            var                           = [P0(i) Ta0(j) condition T0]                                                                                                 ;%[bar,K,-,K] (Another variable that nessery) 
            options                       = odeset('RelTol',1e-3);
            [L ,dep_var]                  = ode45(@func_NiCo,Lspan,dep_var0,options,var) ;

[Walter Roberson]

What happened when you tried with ode45?

Note: in order for us to test we would need your func_NiCo as well

[Jidapa Adam]

When I run ode45 it doesn't stop running and when I change RelTol value to 1e-1 the result is complex number.

load Data_import_NiCo
P0                                    =  [0.5]                                                                                                                        ;%bar         (Intial pressure)
Ta0                                   =  [873]                                                                                                                  ;%K           (Intial air heat temperature)              
repeat                                =  [1]                                                                                                                          ;%Repeat condition Ta constant and Ta non constant
sheet                                 =  0                                                                                                                              ;%Sheet in excel 
condition                             =  0                                                                                                                              ;%Number of condition
for a = 1:length(repeat)  
    for i = 1:length(P0)
        for j = 1:length(Ta0)
            T0                            = Ta0(j)                                                                                                                      ;%K           (Intial temperature before entering bed)
            pj0                           = yj0.*P0(i)                                                                                                                  ;%bar         (Partial pressure)
            Fj0                           = (pj0.*v0_feed)./(R_feed.*T0)                                                                                                ;%mol/s       (Molar flowrate)
            Lspan                         = [0:0.01:0.03]                                                                                                              ;%m           (Bed length range)
            dep_var0                      = [Fj0' T0 Ta0(j)]                                                                                                            ;%[mol/s,K,K] (Intial dependent variable)
            condition                     = condition + 1                                                                                                               ;
            var                           = [P0(i) Ta0(j) condition T0]                                                                                                 ;%[bar,K,-,K] (Another variable that nessery) 
            options                       = odeset('RelTol',1e-3);
            [L ,dep_var]                  = ode45(@func_NiCo,Lspan,dep_var0,options,var) ;

this is func_NiCo

function dydx = func_NiCo(y,x,op)
load Data_import_NiCo
FCH4                    = x(1,1)                                                                                                                    ;
FCO2                    = x(2,1)                                                                                                                    ;
FCO                     = x(3,1)                                                                                                                    ;
FH2                     = x(4,1)                                                                                                                    ;
FH2O                    = x(5,1)                                                                                                                    ;
FC                      = x(6,1)                                                                                                                    ;
FN2                     = x(7,1)                                                                                                                    ;
T                       = x(8,1)                                                                                                                    ;
Ta                      = x(9,1)                                                                                                                    ;
                                                                                                                                                   
P0                      = op(1)                                                                                                                     ;%bar         (Intial pressure)                   
Tac                     = op(2)                                                                                                                     ;%K           (Intial air heat temperature constant)
condition               = op(3)                                                                                                                     ; 
T0                      = op(4)                                                                                                                     ;%K           (Intial temperature)
FT                      = FCH4+FCO2+FCO+FH2+FH2O+FC+FN2                                                                                             ;%mol/s       (Total molar flowrate)                  
%% Partial pressure (bar)
pCH4                    = P0*(FCH4/FT)*(T0/T)                                                                                                       ;
pCO2                    = P0*(FCO2/FT)*(T0/T)                                                                                                       ;
pCO                     = P0*(FCO/FT)*(T0/T)                                                                                                        ;
pH2                     = P0*(FH2/FT)*(T0/T)                                                                                                        ;
pH2O                    = P0*(FH2O/FT)*(T0/T)                                                                                                       ;
pC                      = P0*(FC/FT)*(T0/T)                                                                                                         ;
pN2                     = P0*(FN2/FT)*(T0/T)                                                                                                        ;
%% Heat capacity (J/(mol*K))                                                            
cpj                     = heatcp(:,1) + heatcp(:,2)*T + heatcp(:,3)*(T^2) + heatcp(:,4)*(T^-2)                                                      ;                                                      
%% Del_Heat of reaction
delH1                   = hrxi_Tr(1)+((2.*cpj(3)+2*cpj(4)-cpj(1)-cpj(2)).*(T-Tr))                                                                   ;
delH2                   = hrxi_Tr(2)+((cpj(3)+cpj(5)-cpj(2)-cpj(4)).*(T-Tr))                                                                        ;
delH3                   = hrxi_Tr(3)+((cpj(6)+2.*cpj(4)-cpj(1)).*(T-Tr))                                                                            ;
delH4                   = hrxi_Tr(4)+((cpj(3)+cpj(4)-cpj(5)-cpj(6)).*(T-Tr))                                                                        ;
delH5                   = hrxi_Tr(5)+((2*cpj(3)-cpj(6)-cpj(2)).*(T-Tr))                                                                             ;
%% Reaction rate constant (mol/(kg*s))
ki                      = k0i.*exp(-Eai./T)                                                                                                         ;
%% Equilibrium constant
Kpi                     = Kp0i.*exp(-dHrxi_Tr./T)                                                                                                   ;%[bar^2,-,bar,bar,bar]
%% Adsorption equilibrium constant 
KCH4_1                  = K0j1(1)*exp(hdx1(1)/T)                                                                                                    ;%mol/(kg*s*bar)
KCO2_1                  = K0j1(2)*exp(hdx1(2)/T)                                                                                                    ;%mol/(kg*s*bar)

KCO2_2                  = K0j2(2)*exp(hdx2(2)/T)                                                                                                    ;%bar^-1
KH2_2                   = K0j2(4)*exp(hdx2(4)/T)                                                                                                    ;%bar^-1

KCH4_3                  = K0j3(1)*exp(hdx3(1)/T)                                                                                                    ;%bar^-1
KH2_3                   = K0j3(4)*exp(hdx3(4)/T)                                                                                                    ;%bar^1.5

KCH4_4                  = K0j4(1)*exp(hdx4(1)/T)                                                                                                    ;%bar^-1
KH2_4                   = K0j4(4)*exp(hdx4(4)/T)                                                                                                    ;%bar^1.5
KH2O_4                  = K0j4(5)*exp(hdx4(5)/T)                                                                                                    ;%-

KCO2_5                  = K0j5(2)*exp(hdx5(2)/T)                                                                                                    ;%bar
KCO_5                   = K0j5(3)*exp(hdx5(3)/T)                                                                                                    ;%bar^-1
%% Rate reaction (mol/(kg*s)
r1                      = ((ki(1)*pCH4*pCO2)/(KCO2_1*pCO2+KCH4_1*pCH4))*(1-(((pCO*pH2)^2)/(Kpi(1)*pCH4*pCO2)))                                      ;
r2                      = ((ki(2)*KCO2_2*KH2_2*pCO2*pH2)/((1+KCO2_2*pCO2+KH2_2*pH2)^2))*(1-((pCO*pH2O)/(Kpi(2)*pCO2*pH2)))                          ;   
r3                      = ((ki(3)*KCH4_3)*(pCH4-((pH2^2)/Kpi(3))))/((1+KCH4_3*pCH4+(pH2^1.5/KH2_3))^2)                                              ;
r4                      = ((ki(4)/KH2O_4)*((pH2O/pH2)-(pCO/Kpi(4))))/((1+KCH4_4*pCH4+(pH2O/(KH2O_4*pH2))+(pH2^1.5/KH2_4))^2)                        ;
r5                      = ((ki(5)/(KCO_5*KCO2_5))*((pCO2/pCO)-(pCO/Kpi(5))))/((1+KCO_5*pCO+(pCO2/(KCO_5*KCO2_5*pCO)))^2)                            ;
%% rnet of species (mol/(kg*s))
rn_CH4                  =  -r1-r3                                                                                                                   ; 
rn_CO2                  =  -r1-r2-r5                                                                                                                ;
rn_CO                   =  2*r1+r2+r4+2*r5                                                                                                          ;
rn_H2                   =  2*r1-r2+2*r3+r4                                                                                                          ;
rn_H2O                  =  r2-r4                                                                                                                    ; 
rn_C                    =  r3-r4-r5                                                                                                                 ;
rn_N2                   =  0                                                                                                                        ;
%% 
sumFjcpj                =  FCH4*cpj(1)+ FCO2*cpj(2)+ FCO*cpj(3)+ FH2*cpj(4)+ FH2O*cpj(5)+ FC*cpj(6)+ FN2*cpj(7)                                     ;%Watt/K
sumrihi                 =  r1*-delH1+ r2*-delH2+ r3*-delH3+ r4*-delH4+ r5*-delH5                                                                         ;%Watt/kg   
%% ODE equations 
dFCH4dL                 = Rho_bed*Ac*rn_CH4                                                                                                         ;%mol/(m*s)
dFCO2dL                 = Rho_bed*Ac*rn_CO2                                                                                                         ;%mol/(m*s)
dFCOdL                  = Rho_bed*Ac*rn_CO                                                                                                          ;%mol/(m*s)
dFH2dL                  = Rho_bed*Ac*rn_H2                                                                                                          ;%mol/(m*s)
dFH2OdL                 = Rho_bed*Ac*rn_H2O                                                                                                         ;%mol/(m*s)
dFCdL                   = 0                                                                                                         ;%mol/(m*s)
dFN2dL                  = Rho_bed*Ac*rn_N2                                                                                                          ;%mol/(m*s)      
    if condition <= 6
        dTdL                    = ((Rho_bed*Ac*sumrihi)-(U*ain*Ac*(T-Tac)))/sumFjcpj                                                                ;%K/m       
        dTadL                   = 0                                                                                                                 ;%K/m                    
    else                                                                                                                                                                
        dTdL                    = ((Rho_bed*Ac*sumrihi)-(U*ain*Ac*(T-Ta)))/sumFjcpj                                                                 ;%K/m  
        dTadL                   = (U*aout*Ac_air*(T-Ta))/(mc_air*cpj(8))                                                                            ;%K/m
    end
                                                               
%% Calculations
dydx(1,1)               = dFCH4dL                                                                                                                   ;
dydx(2,1)               = dFCO2dL                                                                                                                   ;
dydx(3,1)               = dFCOdL                                                                                                                    ;
dydx(4,1)               = dFH2dL                                                                                                                    ;
dydx(5,1)               = dFH2OdL                                                                                                                   ;
dydx(6,1)               = dFCdL                                                                                                                     ;
dydx(7,1)               = dFN2dL                                                                                                                    ;
dydx(8,1)               = dTdL                                                                                                                      ;
dydx(9,1)               = dTadL                                                                                                                     ;

end

Answer 1

[Walter Roberson]

You have a "stiff" system of equations. The equations wobble a lot and ode45 has to take very small steps to try to follow the wobbles because they might be important.

[Jidapa Adam]

How can I fix this problem?

[Torsten]

That's what stiff integrators like ODE15S are for - they should be used instead of non-stiff integrators like ODE45 in cases like yours.

[Walter Roberson]

You can fix the problem by stopping interrupting ode45 and let it run the several years it needs to calculate the system properly.

(It is not a bug.)

[Steven Lord]

If you're not familiar with stiffness in ODEs I recommend reading through this post on Cleve Moler's blog. Run at least the one example Cleve recommends running and read his explanation for its performance.