% Modified version of the MATLAB template file SolvIVDifI.m used in the
% solution of Example 30.2 of "A First Course on Kinetics and Reaction
% Engineering."
%
function [nAout,nBout,nYout,nZout,Tout] = Example_30_2_pfr(Tin)
    % Known quantities and constants (in mol m^3 s atm cal K)
    V = 3.0; % m^3
    ntot_in = 1.25; % mol/s
    yAin = 0.5;
    yBin = 0.5;
    P = 2.5; % atm
    cPi = 25.8; % cal/mol/K
    dH = -9800; % cal/mol
    k0 = 2.83e4; % /s
    E = 12900.0; % cal/mol
    R_energy = 1.987; % cal/mol/K
    R_pv = 0.08206 * 1e-3; % m^3*atm/mol/K
    nAin = yAin*ntot_in;
    nBin = yBin*ntot_in;
    nYin = 0;
    nZin = 0;

	% Function that evaluates the ODEs
    function dzdt = odeqns(t,z)
        % t is the value of the independent variable
        % z is a vector containing the values of the dependent variables
        nA = z(1);
        nB = z(2);
        nY = z(3);
        nZ = z(4);
        T = z(5);
        VFR = (nA+nB+nY+nZ)*R_pv*T/P;
        CA = nA/VFR;
        r = k0*exp(-E/R_energy/T)*CA;
        dzdt = [
            -r
            -r
            r
            r
            -dH*r/(nA+nB+nY+nZ)/cPi
        ];
    end % of internal function odeqns
        
    % Initial and final values
    t0 = 0;
    z0 = [
        nAin
        nBin
        nYin
        nZin
        Tin
    ];
    tf = V;
    
    % Solve the ODEs
    [t,zz] = ode45(@odeqns,[t0 tf],z0);

    % Extract the final values of the independent and dependent variables
	last_value = length(t);
	t_f = t(last_value); % final value of the independent variable
	z = zz(last_value,:); % final values of the dependent variables
	
    % set the values of the return variables
    nAout = z(1);
    nBout = z(2);
    nYout = z(3);
    nZout = z(4);
    Tout = z(5);
end % of Example_30_2_pfr.m
