% Modified version of the MATLAB template file SolvNonDif.m used in the
% solution of Example 22.2 of "A First Course on Kinetics and Reaction
% Engineering."
%
function Example_22_2
    % Known quantities and constants (in consistent units)
    CA0 = 10; % mol/gal
    CB0 = 12; % mol/gal
    T0 = 350; % K
    V = 25; % gal
    VFR = 12.5; % gal/min
    dH1_298 = -12; % kJ/mol
    dH2_298 = -21.3; % kJ/mol
    cPA = 85/1000; % kJ/mol/K
    cPB = 125/1000; % kJ/mol/K
    cPD = 200/1000; % kJ/mol/K
    cPU = 170/1000; % kJ/mol/K
    k01 = 1.12e2; % /min
    k02 = 1.87e2; % /min
    E1 = 15.3; % kJ/mol
    E2 = 23.7; % kJ/mol
    R = 8.31446/1000; % kJ/mol/K
    nA0 = CA0*VFR; % mol/min
    nB0 = CB0*VFR; % mol/min
    
    % Function that evaluates the equations
	function f = evalEqns(z)
        % mapping of unknowns into vector z
        nA = z(1);
        nB = z(2);
        nD = z(3);
        nU = z(4);
        T = z(5);
        % calculate variables quantities
        CA = nA/VFR;
        CB = nB/VFR;
        r1 = k01*exp(-E1/R/T)*CA;
        r2 = k02*exp(-E2/R/T)*CB;
        dH1 = dH1_298 + (cPD - cPA - cPB)*(T-298);
        dH2 = dH2_298 + (cPU - cPA - cPB)*(T-298);
        % evaluate the functions being solved
		f = [
			nA0 - nA + V*(-r1 - r2)
            nB0 - nB + V*(-r1 - r2)
            -nD + V*r1
            -nU + V*r2
            (nA0*cPA + nB0*cPB)*(T - T0) + V*(r1*dH1 + r2*dH2)
		];
    end % of internal function evalEqns
    
    % guesses for the solution
    z_guess = [
        nA0
        nB0
        0
        0
        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(5)
    pct_conv_A = 100*(nA0-z(1))/nA0
    sel_D_per_U = z(3)/z(4)
    	
end % of Example_22_2.m