% Modified version of the AFCoKaRE MATLAB template file SolvIVDifD.m used
% in the solution of Activity 28.1 of "A First Course in Kinetics and
% Reaction Engineering."
%
function [conv,vol] = Activity_28_1_pfr
    % Known quantities and constants
    k0 = 1E8; % /min
    E = 50000; % J/mol
    dH = -180000; % J/mol
    cP = 3950; % J/l/K
    T0 = 298; % K
    CA0 = 1.0; % mol/L
    conv = linspace(.2,.9,100)';
    nZ0 = 0; % mol/min
    VFR = 300000/24/60; % L/min
    % Universal constants
    R = 8.3145; % J/mol/K
    % Calculated constant
    nA0 = CA0*VFR;
    
    % Initial values
    t0 = 0;
    z0 = [
        nA0;
        nZ0;
        T0;
    ];
    tf = 1E10;
    
    % loop through the desired conversions
    vol = zeros(100,1);
    fA = 0;
    for i = 1:100
        fA = conv(i);
        % solve the equations
        options = odeset('Events',@stop);
        [t, zz, te, ze, ie] = ode45(@odeqns,[t0, tf],z0,options);
        last_value = length(t);
        vol(i) = t(last_value);
    end
	
	% 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
            -r*dH/VFR/cP;
        ];
    end % of internal function odeqns

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

end % of Activity_28_1_pfr.m
