% Modified version of MATLAB template file SolvIVDifI used to model the PFR
% zones in Activity 35.1
%
function [nAout,nBout,nZout] = Activity_35_pfr(nAin,nBin,nIin,nZin,L)
    % Known quantities and constants (in consistent units)
	D = 0.07; % m
    T = 450 + 273.15; % K
    P = 5.0; % atm
    k = 72.0; % mol/h/atm^0.5/m^3
    R = 8.206E-5; % m^3*atm/mol/K
    
	% Function that evaluates the ODEs
    function dzdt = odeqns(t,z)
        nA = z(1);
        nB = z(2);
        nZ = z(3);
        ntot = nA + nB + nZ + nIin;
        CA = nA*P/ntot/R/T;
        CB = nB*P/ntot/R/T;
        r = k*sqrt(CA)*CB;
        dzdt = [
            -pi()*D^2/2*r
            -pi()*D^2/4*r
            pi()*D^2/2*r
        ];
    end % of internal function odeqns
        
    % Initial and final values
    t0 = 0.;
    z0 = [
        nAin
        nBin
        nZin
    ];
    tf = L;
    
    % Solve the ODEs
    [t,zz] = ode45(@odeqns,[t0 tf],z0);

    % Set the return values
	last_value = length(t);
    nAout = zz(last_value,1);
    nBout = zz(last_value,2);
    nZout = zz(last_value,3);
    
end % of Activity_35_pfr.m
