% Modified version of the AFCoKaRE MATLAB template file SolvIVDifD.m used
% for the solution of Example 29.3 of "A First Course on Kinetics and
% Reaction Engineering."
function [nAout,nYout,nZout,Tout] = Example_29_3_pfr(V,VFRin)
    % Known quantities and constants
    CAin = 1.0; % mol/L
    CYin = 0; % mol/L
    CZin = 0; % mol/L
    Tin = 60 + 273.15; % K
    k0 = 2.63E7; % L/mol/min
    E = 62.; % kJ/mol
    dH = -35.; % kJ/mol
    cP = 0.8; % kJ/L/K
    % Gas constant
    R = 8.3145/1000; % kJ/mol/K
    % Calculated constants
    nAin = VFRin*CAin;
    nYin = VFRin*CYin;
    nZin = VFRin*CZin;

	% Function that evaluates the ODEs
    function dzdt = odeqns(t,z)
        nA = z(1);
        nY = z(2);
        nZ = z(3);
        T = z(4);
        CA = nA/VFRin;
        r = k0*exp(-E/R/T)*CA^2;
        dzdt = [
            -2*r
            r
            r
            -r*dH/VFRin/cP;
        ];
    end % of internal function odeqns
        
    % Initial and final values
    t0 = 0;
    z0 = [
        nAin
        nYin
        nZin
        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);
    nYout = z(2);
    nZout = z(3);
    Tout = z(4);
    
    
end % of Example_29_3_pfr.m
