% Modified version of the MATLAB template file SolvNonDif.m used in the
% solution of Example 30.2 of "A First Course on Kinetics and Reaction
% Engineering."
%
function Equation_30_2
    % Known quantities and constants (in consistent units)
    Ta = 300; % K
    ntot_in = 1.25; % mol/s
    yAin = 0.5;
    yBin = 0.5;
    cPi = 25.8; % cal/mol/K
    nAa = yAin*ntot_in;
    nBa = yBin*ntot_in;
    nYa = 0;
    nZa = 0;
    UA = 13.6;
    
    % Function that evaluates the equations
    function f = evalEqns(z)
        Tb = z(1);
        Td = z(2);
        [nAc,nBc,nYc,nZc,Tc] = Example_30_2_pfr(Tb);
        AMTD = 0.5*((Td-Ta) + (Tc-Tb));
        f = [
            (nAa+nBa+nYa+nZa)*cPi*(Tb-Ta) - (nAc+nBc+nYc+nZc)*cPi*(Tc-Td)
            (nAa+nBa+nYa+nZa)*cPi*(Tb-Ta) - UA*AMTD
        ];
    end % of internal function evalEqns
    
    % guesses for the solution
    z_guess = [
        Ta + 50;
        Ta + 100;
    ];

    % 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)
    
    % Use the value of Tb found by solving the heat exchanger design
    % equations to calculate the outlet molar flow of A and then the
    % conversion
    Tb = z(1);
    Td = z(2);
    [nAc,nBc,nYc,nZc,Tc]=Example_30_2_pfr(Tb);
    display('With thermal backmix:')
    Ta
    Tb
    Tc
    Td
    percent_conversion = 100*(nAa-nAc)/nAa
    
	% Case with no backmix, feed to reactor is at Ta
    [nAb,nBb,nYb,nZb,Tb]=Example_30_2_pfr(Ta);
    display('Without thermal backmix:')
    Tb
    percent_conversion = 100*(nAa-nAb)/nAa

end % of Equation_30_2.m