MATLAB: Error in ode15s (line 150) odeargumen​ts(FcnHand​lesUsed, solver_name, ode, tspan, y0, options, varargin); Error in THT2 (line 46) [v,X] = ode15s(@(v​,x)ODEfun(​v,x,params​),tspan,in​itcons);

ode15s error line 150 line46

I receive the following error but I'm not sure how to go about fixing it. The error message isn't very descriptive:
Error in ode15s (line 150)
    odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
Error in THT2 (line 46)
[v,X] = ode15s(@(v,x)ODEfun(v,x,params),tspan,initcons);
————————————————
This is the code:
function [] = THT2()
clc;
clear all;
close all;
% Define constants
P = 300; % atm
Qf = 0.05713; % m^3 /s
XAf = 0.015;
XNf = 0.985/4;
XHf = 0.985*3/4;
Ac = 1; % m^2
L = 12; % m
Taf = 323; % K

gamma = 0.5; % 1/m
beta = -2.342; % (m^2*s*K) /mol
deltaG = -4.25*10^3; % cal/mol

deltaH = -1.2*10^4; % cal/mol
kminus10 = 7.794*10^11; % atm/s
EoverR = 2*10^4; % K
Rcal = 1.987; % cal /(K*mol)
Ratm = 8.205 * 10^-5; % m^3 atm /(K*mol)
NTf = Qf*P/(Ratm*Taf);     % Total feed molar flux mol/s, N = QP/RT
NAf = NTf*XAf;  % feed molar flux of ammonia
NHf = NTf*XHf;  % feed molar flux of hydrogen
NNf = NTf*XNf;  % feed molar flux of nitrogen
% Length vector
v0 = 0; % meters, top of column, entering the reactor
vf = 12; % meters, bottom of column, extiting reactor core
tspan = linspace(v0, vf, 601);
guess = 1000; % Initial guess of T(z=0).  in K
initcons = [NNf, NHf, NAf, guess, guess]; % initial conditions entering reactor core (z=0)
% 1=N, 2=H, 3=A, 4=T, 5=Ta
params = [P Qf Ac Taf gamma beta deltaG deltaH kminus10 EoverR L Rcal Ratm NTf NAf NHf NNf];
[v,X] = ode15s(@(v,x)ODEfun(v,x,params),tspan,initcons);
subplot(4,1,1);
plot(v, Ta,'r', v, T,'b')
xlabel('Length (m)');
ylabel('Temperature (K)');
title('Temperature versus Length');
legend('Annulus', 'Reactor core');
subplot(4,1,2);
plot(v, NA,'r', v, NN,'b', v, NH,'g')
xlabel('Length (m)');
ylabel('Flux (mol/s)');
title('Concentration vs Length');
legend('Ammonia','Nitrogen','Hydrogen');
subplot(4,1,3);
plot(v, XA)
xlabel('Length (m)');
ylabel('XA');
title('Fractional conversion of ammonia');
subplot(4,1,4);
plot(v, Q)
xlabel('Length (m)')
ylabel ('Flowrate (m^3/s)');
title('Flowrate down the length of the reactor core');
end
function dX = ODEfun (v,x, params);
P       = params(1);
Qf      = params(2);
Ac      = params(3);
Taf     = params(4);
gamma   = params(5);
beta    = params(6);
deltaG  = params(7);
deltaH  = params(8);
kminus10= params(9);
EoverR  = params(10);
L       = params(11);
Rcal    = params(12);
Ratm    = params(13);
NTf     = params(14);
NAf     = params(15);
NHf     = params(16);
NNf     = params(17);
%N2 = x(1),  Molar flowrate of nitrogen
%H2 = x(2),  Molar flowrate of hydrogen
%A = x(3),   Molar flowrate of ammonia
%T = x(4),   Temperature in the reactor core
%Ta = x(5),  Temperature in the annulus
dX = zeros(3,1);
dX(1) = gamma * (x(5) - x(4)); % Temperature profile through annulus = gamma * (Ta - T)
dX(2) = 2*Ac*((kminus10 * exp(EoverR/x(4)))/(Ratm*x(4)*((exp(-deltaG/(Rcal*guess)))^2)*(x(1)*(x(2)^(3/2))/x(3))-(x(3)/(x(2)^(3/2)))));  % Ammonia flux through reactor
dX(3) = -beta*((kminus10 * exp(EoverR/x(4)))/(Ratm*x(4)*((exp(-deltaG/(Rcal*guess)))^2)*(x(1)*(x(2)^(3/2))/x(3))-(x(3)/(x(2)^(3/2))))) + gamma*(x(5) - x(4));  % Temperature profile through reactor core
% krev = (kminus10 * exp(EoverR/x(4)));  % atm/s, solve for k reverse
%K = (exp(-deltaG/(Rcal*guess)));  % equilibrium constant
%K = K*exp(-(deltaH/Rcal)*((1/T)-(1/T))) % From the van't Hoff equation
% Reaction rate equation
%r = (krev/(Ratm*x(4)*(K^2)*(x(1)*(x(2)^(3/2))/x(3))-(x(3)/(x(2)^(3/2)))));
NA = 2*Ac*((kminus10 * exp(EoverR/x(4)))/(Ratm*x(4)*(K^2)*(x(1)*(x(2)^(3/2))/x(3))-(x(3)/(x(2)^(3/2)))));
NN = NNf -.5*(NA-NAf);
NH = NHf - 1.5*(NA-NAf);
XA = NA/(NA+NH+NN); % mole fraction of ammonia leaving the reactor core
Q = (Ratm*x(4)/P)*(NA+NN+NH);
end

Best Answer

It's pretty clear to me:
Undefined function or variable 'guess'.
Error in HTT2>ODEfun (line 103)
    dX(2) = 2*Ac*((kminus10 * exp(EoverR/x(4)))/(Ratm*x(4)*((exp(-deltaG/(Rcal*guess)))^2)*(x(1)*(x(2)^(3/2))/x(3))-(x(3)/(x(2)^(3/2)))));
    % Ammonia flux through reactor
Error in HTT2>@(v,x)ODEfun(v,x,params) (line 44)
    [v,X] = ode15s(@(v,x)ODEfun(v,x,params),tspan,initcons);
Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.
Error in ode15s (line 150)
    odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
Error in HTT2 (line 44)
    [v,X] = ode15s(@(v,x)ODEfun(v,x,params),tspan,initcons);
If you want guess at that point to be the 1000 that you set in your original routine, then you should add it in to how you are handling params . If you want guess to be either the 4th or 5th boundary value, the way you initialize initcons, that would be either x(4) or x(5)
Related Question