% Modified version of the AFCoKaRE MATLAB template file SolvIVDifI.m used
% to model the PFR in Example 31.2 of "A First Course on Kinetics and
% Reaction Engineering."
%
function [nAout,nDout,nUout,Tout,L] = Example_31_2_pfr(nAin,nDin,nUin,Tin)
    % Known quantities and constants (in consistent units)
   k01 = 5.07e6; % /min
   E1 = 65.; % kJ/mol
   k02 = 8.64e4; % dm^3/mol/min
   E2 = 58.1; % kJ/mol
   Cp = 281./1000.; % kJ/dm^3/K
   dH1 = -52.; % kJ/mol
   dH2 = -38.; % kJ/mol
   D = 0.1; % dm
   VFRa = 150.; % dm^3/min
   CAa = 4.0; % mol/dm^3
   RR = 4.0;
   fA = 0.5;
   R = 8.31446/1000.; % kJ/mol/K
   VFRd = VFRa;
   VFRr = RR*VFRd;
   VFRb = VFRr + VFRa;
   nAa = VFRa*CAa;

	% 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);
        nD = z(2);
        nU = z(3);
        T = z(4);
        CA = nA/VFRb;
        k1 = k01*exp(-E1/R/T);
        k2 = k02*exp(-E2/R/T);
        r1 = k1*CA;
        r2 = k2*CA^2;
        dzdt = [
            -pi()*D^2/4*(r1 + r2);
            pi()*D^2/4*r1;
            pi()*D^2/4*r2;
            -pi()*D^2/4*(r1*dH1 + r2*dH2)/VFRb/Cp;
        ];
    end % of internal function odeqns
        
    % Initial and final values
    t0 = 0;
    z0 = [
        nAin
        nDin
        nUin
        Tin
    ];
    
    tf = 1.0E4;
    options = odeset('Events',@stop);
	[t, zz, te, ze, ie] = ode45(@odeqns,[t0, tf],z0,options);

    % Function that provides the integration stopping criterion
    function [stop_when, isterminal, direction] = stop(t,z)
        isterminal = 1;
        direction = 0;
        
        nAd = z(1)/(1+RR);
        stop_when = fA - (nAa - nAd)/nAa;
    end % of internal function stop

    % 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
	
    % Warn if the independent variable reached tf above
    if (t_f >= tf)
        display('WARNING: the final value of the independent variable was')
        display('reached, instead of the final value of the dependent')
        display('variable.');
        display('Set tf to a larger value and re-execute the function')
    end
	
    % Set the return values
    L = t_f;
    nAout = z(1);
    nDout = z(2);
    nUout = z(3);
    Tout = z(4);
    
end % of Example_31_2_pfr.m
