% Modified version of the AFCoKaRE MATLAB Template file SolvNonDif.m used
% in the solution of Example29.2 of "A First Course on Kinetics and
% Reaction Engineering."
%
function [nAout,nDout,nUout,Tout]=Example_29_2_cstr(nAin, nDin, nUin, Tin)
    % Known quantities and constants
    V = 200.0; % L
    VFR = 100.0; % L/min
    dH1 = -21500.0; % cal/mol
    dH2 = -24000.0; % cal/mol
    cP = 1.0*1000; % cal/L/K
    k01 = 3.4E5; % /min
    k02 = 1.67E14; % /min
    E1 = 9000.0; % cal/mol
    E2 = 23000.0; % cal/mol
    % Gas constant
    R = 1.987; % cal/mol/K
    
    % Function that evaluates the equations
    function f = evalEqns(z)
        nA = z(1);
        nD = z(2);
        nU = z(3);
        T = z(4);
        CA = nA/VFR;
        r1 = k01*exp(-E1/R/T)*CA;
        r2 = k02*exp(-E2/R/T)*CA^2;
        f = [
            nAin - nA - V*(r1 + r2)
            nDin - nD + V*r1
            nUin - nU + V*r2
            VFR*cP*(T - Tin) + V*(r1*dH1 + r2*dH2)
        ];
    end % of internal function evalEqns
    
    % guesses for the solution
    z_guess = [
        0.5*nAin
        nDin + 0.5*nAin
        nUin + 0.01*nAin
        Tin + 10
    ];

    % Solve the set of algebraic equations
    options = optimoptions(@fsolve,'MaxFunEvals',2000);
    z = fsolve(@evalEqns, z_guess, options);
    display('The values of the functions being solved are as follows:');
    f = evalEqns(z)
    nAout = z(1);
    nDout = z(2);
    nUout = z(3);
    Tout = z(4);
    
end % of Example_29_2_cstr.m