% Modified version of the AFCoKaRE MATLAB template file SolvNonDif.m used
% in the solution of the second part of Example 29.1 of "A First Course on
% Kinetics and Reaction Engineering."
%
function Example_29_1_b
    % Known quantities and constants
    k0 = 12000.0; % mol/L/min
    E = 6500.0; % cal/mol
    K = 0.21; % mol/L
    VFR = 9.5; % L/min
    CAa = 1.3; % mol/L
    nZa = 0; % mol/min
    Ta = 75 + 273.15; % K
    dH = 4300.0; % cal/mol
    cP = 0.515*1000.0; % cal/L/K
    fA = 0.99;
    % Gas constant
    R = 1.987; % cal/mol/K
    % Calculated constants
    nAa = VFR*CAa;
    nAd = nAa*(1-fA);
    CAd = nAd/VFR;
    
    % Function that evaluates the equations
    function f = evalEqns(z)
        V = z(1);
        nAb = z(2);
        nAc = z(3);
        nZb = z(4);
        nZc = z(5);
        nZd = z(6);
        Tb = z(7);
        Tc = z(8);
        Td = z(9);
        CAb = nAb/VFR;
        CAc = nAc/VFR;
        ra = k0*exp(-E/R/Tb)*CAb/(K + CAb);
        rb = k0*exp(-E/R/Tc)*CAc/(K + CAc);
        rc = k0*exp(-E/R/Td)*CAd/(K + CAd);
        f = [
            nAa - nAb - V*ra
            nZa - nZb + V*ra
            VFR*cP*(Tb - Ta) + V*ra*dH
            nAb - nAc - V*rb
            nZb - nZc + V*rb
            VFR*cP*(Tc - Tb) + V*rb*dH
            nAc - nAd - V*rc
            nZc - nZd + V*rc
            VFR*cP*(Td - Tc) + V*rc*dH
        ];
    end % of internal function evalEqns
    
    % guesses for the solution
    z_guess = [
        50.
        0.3*nAa
        0.6*nAa
        nAa - 0.3*nAa
        nAa - 0.6*nAa
        nAa - nAd
        Ta + 10
        Ta + 20;
        Ta + 30;
    ];

    % Solve the set of algebraic equations
    z = fsolve(@evalEqns, z_guess);
    display('The solver found the following values for the unknowns:');
    z
    display('The corresponding values of the functions being solved are as follows:');
    f = evalEqns(z)
    V_cstr = z(1)
    V_total = 3*z(1)
end % of Example_29_1_b.m