% Modified version of the MATLAB template file SolvIVDifD.m used to model
% the pfr in Example 30.1 of "A First Course on Kinetics and Reaction
% Engineering."
%
function [V,Tout] = Example_30_1_pfr(Tin)
    % Known quantities and constants (in consistent units)
    VFR = 750.0; % L/min
    CAa = 3.8; % mol/L
    dH = -79.8; % kJ/mol
    cP = 987.0 * 4.184e-3; % kJ/L/K
    fA = 0.8;
    k0 = 3.38e6; % /min
    E = 50.0; % kJ/mol
    R = 8.31446e-3; % kJ/mol/K
    nAin = CAa*VFR;
    nZin = 0;

	% Function that evaluates the ODEs
    function dzdt = odeqns(t,z)
        nA = z(1);
        nZ = z(2);
        T = z(3);
        CA = nA/VFR;
        r = k0*exp(-E/R/T)*CA;
        dzdt = [
            -r
            r
            -dH*r/VFR/cP;
        ];
    end % of internal function odeqns
        
    % Initial values
    t0 = 0;
    z0 = [
        nAin
        nZin
        Tin
    ];
     
    tf = 1E6; % This value should NOT be reached
    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;
        
        stop_when = z(1) - nAin*(1-fA);
    end % of internal function stop

    % Extract the outlet values of the independent and dependent variables
	last_value = length(t);
	t_f = t(last_value);
	z = zz(last_value,:);
	
    % 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 value
    V = t_f;
    Tout = z(3);
end % of Example_30_1_pfr.m
