% Modified version of the MATLAB template file SolvNonDif.m used in the
% solution of Example 24.1 of "A First Course on Kinetics and Reaction
% Engineering."
%
function Example_24_1(nAg,nBg,nYg,nZg,Tg)
    % Known quantities and constants (in consistent units)
    V = 500; % cm3
    VFR = 1; % cm3/s
    CA0 = 0.015; % mol/cm3
    CB0 = CA0; % mol/cm3
    nA0 = VFR*CA0; % mol/s
    nB0 = VFR*CB0; % mol/s
    nY0 = 0; % mol/s
    nZ0 = 0; % mol/s
    T0 = 50 + 273.15; % K
    Cp = 0.35*4.184; % J/g/K
    rho = 0.93; % g/cm3
    dH = -20000; % J/mol
    k0 = 3.24e12; % cm3/mol/s (pre-exponential factor in equation (2))
    E = 25000*4.184; % J/mol (activation energy in equation (2))
    R = 8.31446; % J/mol/K
    
    % Function that evaluates the equations
    function f = evalEqns(z)
        % mapping of unknowns into vector z
        nA = z(1);
        nB = z(2);
        nY = z(3);
        nZ = z(4);
        T = z(5);
        % calculate variables quantities
        CA = nA/VFR; % equation (8)
        CB = nB/VFR; % equation (9)
        r1 = k0*exp(-E/R/T)*CA*CB; % equation (2)
        f = [
            nA0 - nA - V*r1
            nB0 - nB - V*r1
            nY0 - nY + V*r1
            nZ0 - nZ + V*r1
            VFR*rho*Cp*(T-T0) + V*r1*dH
        ];
    end % of internal function evalEqns
    
    % guesses for the solution
    z_guess = [
        nAg
        nBg
        nYg
        nZg
        Tg
    ];

    % 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(5)
    pct_conv = 100*(nA0-z(1))/nA0  % equation (10)
    	
end % of Example_24_1.m