% Modified version of the MATLAB template file SolvNonDif.m used in the
% solution of Example 22.1 of "A First Course on Kinetics and Reaction
% Engineering."
%
function Example_22_1(T0)
    % Known quantities and constants (in consistent units)
    VFR = 0.08; % ft3/min
    nA0 = 0.01; % lbmol/min
    nB0 = 0.25; % lbmol/min
    nZ0 = 0; % lbmol/min
    V = 18; % ft3
    CpA = 1000; % BTU/lbmol/degR
    CpB = 180; % BTU/lbmol/degR
    CpZ = 1200; % BTU/lbmol/degR
    dH = -1.7e4; % BTU/lbmol
    k0 = 1.2e14; % ft3/lbmol/min
    EoverR = 23000; % degR
    K0 = 6.5e-13; % ft3/lbmol
    HoverR = 20000; % degR
    
    % Function that evaluates the equations
    function f = evalEqns(z)
        % mapping of unknowns into vector z
        nA = z(1);
        nB = z(2);
        nZ = z(3);
        T = z(4);
        % calculate variables quantities
        CA = nA/VFR;
        CB = nB/VFR;
        CZ = nZ/VFR;
        r1 = k0*exp(-EoverR/T)*CA*CB*(1-CZ/(K0*exp(HoverR/T)*CA*CB));
        f = [
            nA0 - nA - V*r1
            nB0 - nB - V*r1
            nZ0 - nZ + V*r1
            (nA0*CpA + nB0*CpB + nZ0*CpZ)*(T - T0) + V*r1*dH
        ];
    end % of internal function evalEqns
    
    % guesses for the solution
    z_guess = [
        nA0
        nB0
        nZ0
        T0
    ];

    % Solve the set of algebraic equations and display the results
    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)
  
    % Report the results
    T = z(4)
    pct_conv = 100*(nA0-z(1))/nA0
    	
end % of Example_22_1.m