% Modified version of the AFCoKaRE MATLAB template file SolvIVDifI.m used
% to model the PFR in Activity 31.1 of "A First Course on Kinetics and
% Reaction Engineering."
%
function [nAout,nBout,Tout] = Activity_31_1_pfr(nAin,nBin,Tin,VFRin)
    % Known quantities and constants
    dH = -14.0; % kcal/mol
    k0 = 4.2E15; % cc/mol/min
    E = 18.0; % kcal/mol
    cP = 1.3/1000.; % kcal/cc/K
    D = 5.0; % cm
    L = 50; % cm
    R = 1.987/1000.; % kcal/mol/K
    % Calculated constants
    Axs = pi()*D^2/4.0; % cm^2
    
	% 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);
        T = z(3);
        CA = nA/VFRin;
        CB = nB/VFRin;
        rate = k0*exp(-E/R/T)*CA*CB;
        dzdt = [
            -Axs*rate;
            Axs*rate;
            -Axs*rate*dH/VFRin/cP;
        ];
    end % of internal function odeqns
        
    % Initial and final values
    t0 = 0;
    z0 = [
        nAin;
        nBin;
        Tin
    ];
    tf = L;
    
    % 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 return values
    nAout = z(1);
    nBout = z(2);
    Tout = z(3);
    
end % of Activity_31_1_pfr.m
