% Modified version of the MATLAB template file SolvIVDifI.m used in the
% solution of Example 26.3 of "A First Course on Kinetics and Reaction
% Engineering."
%
function Example_26_3
    % Known quantities and constants
    D = 10; % cm
    L = 500; % cm
    CA0 = 1.0/1000.0; % mol/cm^3
    CB0 = 1.2/1000.0; % mol/cm^3
    nY0 = 0;
    nZ0 = 0;
    T0 = 30 + 273.15; % K
    k0 = 8.72E5*1000.; % cm^3/(mol min)
    E = 7200; % cal/mol
    dH = -10700; % cal/mol
    Cp = 1.0; % cal/(g K)
    rho = 1.0; % g/cm^3
    VFR = 75*1000; % cm^3/min
    % Universal constants
    R = 1.987; % cal/(mol K)
    % Calculated constants
    nA0 = CA0*VFR;
    nB0 = CB0*VFR;
    
    % Function that evaluates the ODEs
    function dzdt = odeqns(t,z)
        nA = z(1);
        nB = z(2);
        nY = z(3);
        nZ = z(4);
        T = z(5);
        CA = nA/VFR; % equation (10)
        CB = nB/VFR; % equation (11)
        r = k0*exp(-E/R/T)*CA*CB; % equation (2)
        dzdt = [
            -pi()*D^2/4*r % equation (4)
            -pi()*D^2/4*r % equation (5)
            pi()*D^2/4*r % equation (6)
            pi()*D^2/4*r % equation (7)
            -pi()*D^2/4*r*dH/VFR/rho/Cp % equation (9)
        ];
    end % of internal function odeqns

    % Initial and final values
    t0 = 0;
    z0 = [
        nA0
        nB0
        nY0
        nZ0
        T0
    ];
    tf = L;

    % 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,:);
    
    % Generate the plots
    figure
    plot(t,zz(:,1))
    title('Molar Flow of A vs. Axial Position')
    xlabel('Axial Position z (m)')
    ylabel('Molar Flow Rate of A (mol/min)')
    figure
    plot(t,zz(:,5)-273.15)
    title('Temperature vs. Axial Position')
    xlabel('Axial Position z (m)')
    ylabel('Temperature (deg C)')
    
end % of Example_26_3.m
