How to frequently run ODE45?

Hi everyone. I have written the following code in order to solve a set of ODE equations. As you can see in my function dy=equation(z,y) I have a constant named T (T=300 %temperatue).

 function bimolecular_reaction
function Kads = adsparameter(P,A,T,m)
 Kb = 1.3806488e-23; %Boltzmann Constant 
(J/K)
 Kads = (P*A)/(sqrt(2*pi*m*Kb*T));
 end
 function Kdes = 
desparameter(A,T,m,tetarot,sigma,Edes)
 Kb = 1.3806488e-23; %Boltzmann Constant 
(J/K)
 h = 6.2606957e-34; %Planck Constant (J/s)
 R = 8.3144621; %universal gas constant 
(J/mol.K)
 Kdes = 
((Kb*(T^3))/(h^3))*((A*2*pi*m*Kb)/(sigma*tetarot))*
exp((-Edes)/(R*T));
 end
 function Karr = surfaceparameter (T,Eac,nu)
 %Kb = 1.3806488e-23; %Boltzmann Constant 
(J/K)
 %h = 6.2606957e-34; %Planck Constant (J/s)
 R = 8.3144621; %universal gas constant 
(J/mol.K)
 Karr = nu*exp((-Eac)/(R*T));
 end
 function dy = equations(z,y)
 T=300; %temperature
 P = 20e5; %total presure of system (Pa)
 pa = (2/3)*P; %partial pressure of CO (Pa)
 pb = (1/3)*P; %partial pressure of O2 (Pa)
 pc = 0; %partial pressure CO2 (pa)
 ma = 28 * 1.66054e-27; % mass of CO (Kg)
 mb = 32 * 1.66054e-27; %mass of O2 (Kg)
 mc = 80 * 1.66054e-27; %mass of CO2 (Kg)

 k1ads = adsparameter(pa,1e-20,T,ma)*1e-13;
k1des = desparameter(1e-20,T,ma,2.8,1,80e3)*1e-
13;
 k2ads = adsparameter(pb,1e-20,T,mb)*1e-13;
 k2des = desparameter(1e-20,T,mb,2.08,2,40e3)*1e-
13;
 kf = surfaceparameter(T,120e3,1e13)*1e-13;
 kb = surfaceparameter(T,180e3,1e13)*1e-13;
 k4ads = adsparameter(pc,1e-20,T,mc)*1e-13;
 k4des = desparameter(1e-
20,T,mc,0.561,1,10e3)*1e-13;
 r1f = k1ads*y(4);
 r1b = k1des*y(3);
 r2f = k2ads*(y(4)^2);
 r2b = k2des*(y(2)^2);
 r3f = kf*y(1)*y(2);
 r3b = kb*y(3)*y(4);
 r4f = k4ads*y(4);
 r4b = k4des*y(3);
 dy = zeros (4,1);
 dy(1) = r1f-r1b-r3f+r3b;
 dy(2) = (2*r2f)-(2*r2b)-r3f+r3b;
 dy(3) = r4f-r4b+r3f-r3b;
 dy(4) = r1b-r1f-(2*r2f)+(2*r2b)+r3f-r3b-r4f+r4b;
 end
[z,y] = ode45(@equations,[0 1],[0 0 0 1]);
 rco2 = desparameter(1e-20,300,80 * 1.66054e-
27,0.561,1,10e3)*y(end,3);
 plot(z,y(:,2),'r')
end

Firstly, I am going to construct a loop by which I could solve my ode system based on different values of T (between 300 to 1000) but I do not know how can do so. Secondly, since, the constants of my ode system is of order of 10^13 it will takes me too much time to reach to an answer for my ode system how can I also solve this problem?

[Walter Roberson]

http://www.mathworks.com/help/matlab/math/parameterizing-functions.html for the matter of how to pass in T values.

I wonder if you should be using one of the "stiff" solvers such as ode23s ?

[Hamid Ahmadi Eshtehardi]

Thank's for your comment. I have already read the parameterzing functions but I did not understand that how can I use them in my code. Could you please give me a little bit more advise in this? About the solver you were right my system was stiff so I had to use one of stiff solvers such as ode15s or ode23s.

Answer 1

[Walter Roberson]

for T = 300:100:1000
   [z,y] = ode23s(@(t,y) equations(t,y,T), [0 1], [0 0 0 1]);
   rco2 = desparameter(1e-20,T,80 * 1.66054e-27,0.561,1,10e3)*y(end,3);
   plot(z,y(:,2),'r', 'DisplayName', sprintf('T=%g', T));
   legend show
   hold on
   drawnow();
end

 function dy = equations(z,y,T)
 %T=300; %temperature
 etc