% Modified version of the MATLAB template file SolvIVDifI.m used to model
% the PFR in Example 33.1 of "A First Course on Kinetics and Reaction
% Engineering."
%
function [PFR_length,PFR_conv] = Example_33_1_pfr
    % Known quantities and constants (in consistent units)
    k = 1.2E-2; % /s
    CAfeed = 1.0e-3; % mol/m^3
    Peax = 2.0;
    us = 0.01; % m/s
    d = 0.004; % m
    L = 50*d;

	% Function that evaluates the ODEs
    function dzdt = odeqns(t,z)
        CA = z(1);
        r = k*CA;
        dzdt = [
            -r/us;
        ];
    end % of internal function odeqns
        
    % Initial and final values
    t0 = 0.0;
    z0 = [
        CAfeed
    ];
    tf = L;
    
    % 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
	
    PFR_length = t;
    PFR_conv = (CAfeed - zz(:,1))*100/CAfeed;
    
end % of Example_33_1_pfr.m
