function dx = SMR_method(z,x)
% variables
xch4 = x(1); xh2o = (2); xh2 = (3); xco = (4); xco2 = (5);
% Density of catalyst (kg/m^3)
y = 2355.2;
% molar feed rate of ch4 (kmole/hr)
Fch4 = 5.1688;
% Temperature (k)
T = 793.15;
% Cross-sectional area (m^2)
A = 0.0324;
% Gas constant J/mole.K
R = 8.314;
%pressure bar
P = 29;
% rate constant k1,k2,k3 : bar^0.5,bar^-1,bar^0.5
k1 = 4.225e15*exp(-240100/R*T);
k2 = 1.955e6*exp(-67130/R*T);
k3 = 1.020e15*exp(-243900/R*T);
% Equilibrium constant K1,K2,K3 : bar^2,bar^0,bar^2
K1 = exp(-26830/(T+30.114));
K2 = exp(4400/(T-4.036))*(-1);
K3 = K1*K2;
% Adsorption constants Kch4,Kh2,kco,kh20 : bar^-1,bar^-1,bar^-1,bar^0
Kch4 = 6.65e-4*exp(38280/R*T);
Kh2 = 6.12e-9*exp(82900/R*T);
Kco = 8.23e-5*exp(70650/R*T);
Kh2o = 1.77e5*exp(-88680/R*T)*(-1);
% reaction rate
r1 = k1(xch4*xh2o*P^1.5/((xh2)^0.5)-(xh2)^2.5*xco*P^3.5/K1)/(1+Kh2o*xh2o/xh2+P(Kch4*xch4+Kco*xco+Kh2*xh2))^2;
r2 = k2(xco*xh2o*P/xh2-xco2/K2)/(1+Kh2o*xh2o/xh2+P(Kch4*xch4+Kco*xco+Kh2*xh2))^2;
r3 = k3(xch4*xh2o^2/(P^0.5*(xh2)^3.5)-(xh2)^0.5*xco2*P^1.5/K3)/(1+Kh2o*xh2o/xh2+P(Kch4*xch4+Kco*xco+Kh2*xh2))^2;
% mass balance
dxch4 = A*y*(-r1-r3)/Fch4;
dxh20 = A*y*(-r1-r2-2*r3)/Fch4;
dxh2 = A*y*(3*r1+r2+4*r3)/Fch4;
dxco = A*y*(r1-r2)/Fch4;
dxco2 = A*y*(r2+r3)/Fch4;
% Assign output variables
dx(1,:) = dxch4;
dx(2,:) = dxh20;
dx(3,:) = dxh2;
dx(4,:) = dxco;
dx(5,:) = dxco2;
end
% Initial mass fraction
x0 = [0.2128, 0.7147, 0.0259, 0, 0.0119];
% reactor length (m)
h = [0, 12];
% Run ODE solver
[z, y] = ode23s(@SMR_method, h, x0);
>> SMR_runfile
Array indices must be positive integers or logical values.
Error in SMR_method (line 49)
r1 = k1(xch4*xh2o*P^1.5/((xh2)^0.5)-(xh2)^2.5*xco*P^3.5/K1)/(1+Kh2o*xh2o/xh2+P(Kch4*xch4+Kco*xco+Kh2*xh2))^2;
Error in odearguments (line 92)
f0 = ode(t0,y0,args{:}); % ODE15I sets args{1} to yp0.
Error in ode23s (line 122)
= odearguments(odeIsFuncHandle, odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);
Error in SMR_runfile (line 14)
[z, y] = ode23s(@SMR_method, h, x0);
please help
