MATLAB: Using for loop with ODE solver

changing parametersode's

Hello MATLAB central!
I have written a function that reads Arrhenius equation and uses different parameters in different temperature ranges, like so:
 function dxdT_MSK =MSK_pyrolysis(T_MSK,x_MSK)
    for k = 1:421
    if T_MSK(k,1) < 500
        dxdT_MSK(k,1) = A_MSK(1).*exp(-E_MSK(1)./(R.*T_MSK(k,1))).*(1-x_MSK(k,1))^n_MSK(1);
    elseif (T_MSK(k,1) >= 500) && (T_MSK(k,1) <620)
               dxdT_MSK(k,1) = A_MSK(2).*exp(-E_MSK(2)./(R.*T_MSK(k,1))).*(1-x_MSK(k,1))^n_MSK(2);
    elseif (T_MSK(k,1) >= 620) && (T_MSK(k,1) <690)
        dxdT_MSK(k,1)=A_MSK(3).*exp(-E_MSK(3)./(R.*T_MSK(k,1))).*(1-x_MSK(k,1))^n_MSK(3); 
    elseif T_MSK(k,1) >=690
        dxdT_MSK(k,1)=0;
    end
    end
 end
I would like to solve it with an ODE solver. I have attempted to do it by using this command:
 T_span_MSK = linspace(460,860)';
 E_MSK = []; %matrix of acitvation energies for pure MSK (primary, secondary, and teriary degradations respectively)
n_MSK = []; %matrix of reaction orders for pure MSK
A_MSK =[]; %matrix of preexponential factors for pure MSK
E_MSK(1) = 8.309;
n_MSK(1) =3;
A_MSK(1) = 9.76*10^(-2);
E_MSK(2) = 37.832;
n_MSK(2) = 0.9;
A_MSK(2) = 3.64*10^2;
E_MSK(3)=14.99;
n_MSK(3) = 1;
A_MSK(3) = 1.61;
[T_MSK,x_MSK] = ode45(@MSK_pyrolysis, T_span_MSK,0);
When I attempt to do so, I get this error message:
Index in position 1 exceeds array bounds (must not exceed 1).
Error in MSK_pyrolysis (line 6)
    if T_MSK(k,1) < 500
Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.
Error in ode45 (line 115)
  odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
Error in Untitled2 (line 13)
[T_MSK,x_MSK] = ode45(@MSK_pyrolysis, T_span_MSK,0);
I would really appreciate some pointers as to what I could improve or change in my code. I know it would be possible to just do it seperately for each temperature range and then use "brute force" to combine them together, but obviously there must be some more sleek way of doing it.
KInd regards.

Best Answer

Structurally, your code needs to look more like the following:
 T_span_MSK = linspace(460,860)';
%  E_MSK = []; %matrix of acitvation energies for pure MSK (primary, secondary, and teriary degradations respectively)
% n_MSK = []; %matrix of reaction orders for pure MSK
% A_MSK =[]; %matrix of preexponential factors for pure MSK
E_MSK(1) = 8.309;
n_MSK(1) =3;
A_MSK(1) = 9.76*10^(-2);
E_MSK(2) = 37.832;
n_MSK(2) = 0.9;
A_MSK(2) = 3.64*10^2;
E_MSK(3)=14.99;
n_MSK(3) = 1;
A_MSK(3) = 1.61;
[T_MSK,x_MSK] = ode45(@(T_MSK,x_MSK) MSK_pyrolysis(T_MSK,x_MSK,A_MSK,E_MSK,n_MSK), T_span_MSK,0);
 function dxdT_MSK =MSK_pyrolysis(T_MSK,x_MSK, A_MSK,E_MSK,n_MSK)
    R = 8.3145; % J mol^-1 K^-1
    if T_MSK < 500
        dxdT_MSK = A_MSK(1).*exp(-E_MSK(1)./(R.*T_MSK)).*(1-x_MSK)^n_MSK(1);
    elseif (T_MSK >= 500) && (T_MSK <620)
               dxdT_MSK = A_MSK(2).*exp(-E_MSK(2)./(R.*T_MSK)).*(1-x_MSK)^n_MSK(2);
    elseif (T_MSK >= 620) && (T_MSK <690)
        
        dxdT_MSK=A_MSK(3).*exp(-E_MSK(3)./(R.*T_MSK)).*(1-x_MSK)^n_MSK(3); 
    elseif T_MSK >=690
        dxdT_MSK=0;
    end
   
 end
However,
(i) when x_MSK gets to 1 you start getting complex values because of the value of n_MSK(2).
(ii) you don't give units so I might have used the wrong value for R (which you didn't include).
(iii) you don't have a beta in your rate equations (dx.dT = A/beta*exp(-E/(RT))*(1-x)^n)- is that because that constant is wrapped up in your A values?
Related Question