% Modified version of the AFCoKaRE MATLAB template file SolvNonDif.m used
% to model the recycle PFR in Example 31.2 of "A First Course on Kinetics 
% and Reaction Engineering."
%
function Example_31_2
    % Known quantities and constants (in consistent units)
   Cp = 281./1000.; % kJ/dm^3/K
   VFRa = 150.; % dm^3/min
   Ta = 420; % K
   CAa = 4.0; % mol/dm^3
   RR = 4.0;
   fA = 0.5;
   R = 8.31446/1000.; % kJ/mol/K
   VFRd = VFRa;
   VFRr = RR*VFRd;
   VFRb = VFRr + VFRa;
   VFRc = VFRb;
   nAa = VFRa*CAa;
   nDa = 0;
   nUa = 0;

    % Function that evaluates the equations
    function f = evalEqns(z)
        nAb = z(1);
        nDb = z(2);
        nUb = z(3);
        Tb = z(4);
        [nAc,nDc,nUc,Tc,L] = Example_31_2_pfr(nAb,nDb,nUb,Tb);
        f = [
            nAa + RR*nAc/(1+RR) - nAb;
            nDa + RR*nDc/(1+RR) - nDb;
            nUa + RR*nUc/(1+RR) - nUb;
            VFRa*Cp*(Tb-Ta) + VFRr*Cp*(Tb-Tc);
        ];
    end % of internal function evalEqns
    
    % guesses for the solution
    z_guess = [
        0.5*nAa*(1+RR);
        0.1*nAa*(1+RR);
        0.1*nAa*(1+RR);
        Ta + 20;
    ];

    % 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)
    
    % Solve the PFR design equations
    nAb = z(1);
    nDb = z(2);
    nUb = z(3);
    Tb = z(4);
    [nAc,nDc,nUc,Tc,L] = Example_31_2_pfr(nAb,nDb,nUb,Tb);
    nDd = nDc/(1+RR);
    nUd = nUc/(1+RR);
    L
    sel_D_per_U = nDd/nUd
    Td = Tc
	
end % of Example_31_2.m