% Modified version of the AFCoKaRE MATLAB Template file SolvIVDifI.m used
% in the solution of Example 29.2 of "A First Course on Kinetics and
% Reaction Engineering."
%
function [nAout,nDout,nUout,Tout] = Example_29_2_pfr(nAin,nDin,nUin,Tin)
    % Known quantities and constants
    V = 200.0; % L
    VFR = 100.0; % L/min
    dH1 = -21500.0; % cal/mol
    dH2 = -24000.0; % cal/mol
    cP = 1.0*1000; % cal/L/K
    k01 = 3.4E5; % /min
    k02 = 1.67E14; % /min
    E1 = 9000.0; % cal/mol
    E2 = 23000.0; % cal/mol
    % Gas constant
    R = 1.987; % cal/mol/K

	% Function that evaluates the ODEs
    function dzdt = odeqns(t,z)
        nA = z(1);
        nD = z(2);
        nU = z(3);
        T = z(4);
        CA = nA/VFR;
        r1 = k01*exp(-E1/R/T)*CA;
        r2 = k02*exp(-E2/R/T)*CA^2;
        dzdt = [
            -r1 - r2
            r1
            r2
            -(r1*dH1 + r2*dH2)/VFR/cP
        ];
    end % of internal function odeqns
        
    % Initial and final values
    t0 = 0;
    z0 = [
        nAin
        nDin
        nUin
        Tin
    ];
    tf = V;
    
    % Solve the ODEs
    [t,zz] = ode45(@odeqns,[t0 tf],z0);

    % Return the values of the independent and dependent variables when the
    % integration stopping condition was satisfied
	last_value = length(t);
	t_f = t(last_value);
	z = zz(last_value,:);
    nAout = z(1);
    nDout = z(2);
    nUout = z(3);
    Tout = z(4);
	
end % of Example_29_2_pfr.m
