MATLAB: Way to visualise system of ODEs

dsolveMATLABmodelode45

Hi,
I have a problem involving a Plug-Flow reactor where I have a system of coupled differential equations both with respect to the length travelled through the reactor.
A link to the theory; Here
I essetially have two ODEs that need to be solved;
dX/dL = n*R*A/(2*Fn)
dT/dL = n*H*A*R/(Ft*Cp)
where; H=f(T), R = f(T), Cp = f(T), n = f(T,X) etc. My code is included below.
I believe I cannot solve this analytically (dsolve) so was thinking down the lines of ode45 etc, however I'm not sure why my results so far haven't been successful.
Any help with where to begin would be amazing!
    clear; close all; format shortg; opengl software; clc; 
    %% Input specifications
    R = 8.314;                                                                 %Gas Const [J/molK]
    alpha = 0.5;                                                               %Constant for use in Rate Eqn
    D = 1;                                                                     %Diameter of Reactor
    A = pi*D^2/4;                                                              %Area of Reactor
    % Temp & Pressure
    T_inC = 306;                                                               %[DegC]
    T_inK = T_inC+273.15;                                                      %[K]
    P_in = 208;                                                                %[bar] ----- subject to change
    P_inA = P_in*1.01325;                                                      %[atm] -------||--------------
    % Initial Mole Fractions 
    F_H2_in = 73.95; F_N2_in = 20.86; F_NH3_in = 1.24; 
    F_T_in = F_N2_in+F_H2_in+F_NH3_in;
    
    Y_N2 = F_N2_in/F_T_in; Y_H2 = F_H2_in/F_T_in; 
    Y_NH3 = F_NH3_in/F_T_in;
    sumY = Y_N2+Y_H2+Y_NH3;                                                    %Make sure sum(Yi)=1
    %% Equilibrium
    disp('------------------------------------------')
    disp('Starting Simulation');
    syms T(L) X(L)
    
    P = P_inA;
    log10Ka = -2.691122.*log(T)-5.519265e-5.*T...
        +1.848863e-6.*T.^2+2001.6./T+2.689; 
    Ka = 10.^(log10Ka);
    % Arrhenius
    A0 = 8.849e14;                                                             %Arrhenius Constant
    Ea = 170560;                                                               %Activation Energy [J/mol]
    k = A0.*exp(Ea./(R.*T));                                                   %Rate Constant
    
    % Fugacities & Activities
    phi_N2 = 0.934+0.203e-3.*T+0.296e-3*P-...
        0.271e-6.*T.^2+0.4775e-6*P^2;                                          %Fugacity Coeff - N2 (Reactant)
    y_N2 = Y_N2*(1-X)./(1-2.*X.*Y_N2);                                         %Mol of N2
    a_N2 = phi_N2.*y_N2.*P;                                                    %Activity of N2


    phi_H2 = exp(exp(-3.84.*T.^0.125+0.541)*P...
        -exp(-0.126.*T.^0.5-15.98)*P^2+...
        300*exp(-0.0119.*T-5.941)*exp(+P/300));                                %Fugacity Coeff - H2 (Reactant)
    y_H2 = (Y_H2-3.*X.*Y_N2)./(1-2.*X.*Y_N2);                                  %Mol of H2
    a_H2 = phi_H2.*y_H2.*P;                                                    %Activity of N2
    phi_NH3 = 0.144+0.203e-2.*T-0.449e-3*P-...
        0.114e-5.*T.^2+0.276e-6*P^2;                                           %Fugactiy Coeff - NH3 (Product)
    y_NH3 = (Y_NH3+2.*X.*Y_N2)./(1-2.*X.*Y_N2);                                %Mol of NH3
    a_NH3 = phi_NH3.*y_NH3.*P;                                                 %Activity of N2
    % Other Vars
    C_pT = 35.31+0.02*T+0.00000694*T^2-0.0056*P+0.000014*P*T;                  %Specific Heat Cap. of Mixture
    F_N2 = F_N2_in*(1-X); 
    F_H2 = F_H2_in-3*X*F_N2_in;
    F_NH3 = F_NH3_in+2*X*F_N2_in;
    F_T = F_N2+F_H2+F_NH3;
    
    % Final Vars
    dH_R = 4.184*(-(0.54526+840.609./T+459.734e7./(T.^3))*P-...
        5.34685.*T-0.2525e-3.*T.^2+1.69167e-6.*T.^3-9157.09);                  %Heat of Reaction
    R_NH3 = 2.*k.*(Ka.^2.*a_N2.*(a_H2.^3./a_NH3.^2).^alpha-...
        (a_NH3.^2./a_H2.^3).^(1-alpha));                                       %Rate of Reaction
    eps = -17.539096+0.07697849.*T+6.900548.*T.^2-1.082e-4.*T.^3-...
        26.4247.*T.^4+4.9276e-8.*T.^5+38.937.*X.^3;                            %Catalyst Effectiveness
    
    disp('Variables Simulated');
    
    
    %% Solve Differential Equations
    disp('Starting ODE solving');
    ODE1 = diff(X) == eps.*R_NH3.*A./(2.*F_N2);
    ODE2 = diff(T) == eps.*(-dH_R).*A.*R_NH3./(F_T.*C_pT);
    ODEs = [ODE1; ODE2];
    V1 = odeToVectorField(ODE1,ODE2); 
    F1 = matlabFunction(V1,'vars',{'L','Y','X','T'});
    cond1 = T(0) == T_inK;
    cond2 = X(0) == 0;
    conds = [cond1; cond2];
    Lspan = [0, 30];
    y0 = [T_inK, 0];
    s = ode45(F1)
    disp('ODEs Solved');

Best Answer

Try your code with these changes: 
    [V1,Sbs] = odeToVectorField(ODE1,ODE2); 
    F1 = matlabFunction(V1,'vars',{'t','Y'});
    cond1 = T(0) == T_inK;
    cond2 = X(0) == 0;
    conds = [cond1; cond2];
    Lspan = [0, 30];
    y0 = [T_inK, 0];
    [t,ys] = ode45(F1, Lspan, y0);
    disp('ODEs Solved');
    
figure
yyaxis left
plot(t,ys(:,1))
yyaxis right
plot(t, ys(:,2))
legend(string(Sbs), 'Location','E')
grid
The code prior to that is unchanged. The first argument to ‘F1’ must be the independent variable vector, even if you don’t use it. Since the only variable in ‘F1’ is ‘Y’ otherwise (as the default dependent variable), that’s all matlabFunction needs to know.
Related Question