% Modified version of the AFCoKaRE MATLAB template file SolvNonDif.m used
% to model the recycle PFR in Activity 31.1 of "A First Course on Kinetics 
% and Reaction Engineering."
%
function Activity_31_1
    % Known quantities and constants (in consistent units)
    CAa = 2.0/1000.; % mol/cc
    VFRa = 500.; % cc/min
    Ta = 300.; %K
    RR = 1.3;
    cP = 1.3/1000.; % kcal/cc/K
    % Calculated constants
    nAa = CAa*VFRa;
    nBa = 0;
    VFRd = VFRa;
    VFRr = RR*VFRd;
    VFRb = VFRr + VFRa;
    VFRc = VFRb;

    % Function that evaluates the equations
    function f = evalEqns(z)
        nAb = z(1);
        nBb = z(2);
        Tb = z(3);
        [nAc,nBc,Tc] = Activity_31_1_pfr(nAb,nBb,Tb,VFRb);
        f = [
            nAa + RR*nAc/(1+RR) - nAb;
            nBa + RR*nBc/(1+RR) - nBb;
            VFRa*cP*(Tb-Ta) + VFRr*cP*(Tb-Tc)
        ];
    end % of internal function evalEqns
    
    % guesses for the solution
    z_guess = [
        0.9*CAa*VFRa*(1+RR)
        0.1*CAa*VFRa*(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)
    nBb = z(2)
    Tb = z(3)
    [nAc,nBc,Tc] = Activity_31_1_pfr(nAb,nBb,Tb,VFRb)
    nAd = nAc/(1+RR)
    nBd = nBc/(1+RR)
    Td = Tc
    CAd = nAd/VFRd*1000 % mol/L
    CBd = nBd/VFRd*1000 % mol/L
	
end % of Activity_31_1.m