% Modified version of the AFCoKaRE MATLAB template file SolvIVDifI.m used
% to model the PFR in Example 31.1 of "A First Course on Kinetics and
% Reaction Engineering."
%
function [mSout,mXout] = Example_31_1_pfr(mSin,mXin,VFRin)
    % Known quantities and constants (in consistent units)
    V = 4000.0; % cc
    k1 = 0.034; % /min
    k2 = 0.001; % g/cc

	% 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
        mS = z(1);
        mX = z(2);
        CS = mS/VFRin;
        CX = mX/VFRin;
        rate = k1*CS*CX/(k2 + CS);
        dzdt = [
            -2.2*rate;
            rate;
        ];
    end % of internal function odeqns
        
    % Initial and final values
    t0 = 0;
    z0 = [
        mSin;
        mXin;
    ];
    tf = V;
    
    % 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
    mSout = z(1);
    mXout = z(2);
    
end % of Example_31_1_pfr.m
