% Modified version of the MATLAB template file SolvIVDifI.m used in the
% solution of Activity 23.2 of "A First Course on Kinetics and Reaction
% Engineering."
%
function Activity_23_2
    % Known quantities and constants (in consistent units)
    VFR = 0.16; % ft3/min
    nA0 = 0.02; % lbmol/min
    nB0 = 0.5; % lbmol/min
    nZ0 = 0; % lbmol/min
    T0 = 650; % degR
    V = 18; % ft3
    CpA = 1000; % BTU/lbmol/degR
    CpB = 180; % BTU/lbmol/degR
    CpZ = 1200; % BTU/lbmol/degR
    dH = -1.7e4; % BTU/lbmol
    k0 = 1.2e14; % ft3/lbmol/min
    EoverR = 23000; % degR
    K0 = 6.5e-13; % ft3/lbmol
    HoverR = 20000; % degR
    nAout0 = 0.0005; % lbmol/min
    nBout0 = 0.2405; % lbmol/min
    nZout0 = 0.0095; % lbmol/min
    Tout0 = 653; % degR
    
	% Function that evaluates the ODEs
    function dzdt = odeqns(t,z)
        nA = z(1);
        nB = z(2);
        nZ = z(3);
        T = z(4);
        % calculate variables quantities
        CA = nA/VFR;
        CB = nB/VFR;
        CZ = nZ/VFR;
        r1 = k0*exp(-EoverR/T)*CA*CB*(1-CZ/(K0*exp(HoverR/T)*CA*CB));
        % calculate the sums for the energy balance
        topTerm = (nA0*CpA + nB0*CpB + nZ0*CpZ)*(T - T0);
        bottomTerm = V/VFR*(nA*CpA + nB*CpB + nZ*CpZ);
        dzdt = [
            VFR/V*(nA0 - nA - V*r1)
            VFR/V*(nB0 - nB - V*r1)
            VFR/V*(nZ0 - nZ + V*r1)
            (-topTerm - V*r1*dH)/bottomTerm;
        ];
    end % of internal function odeqns
        
    % Initial and final values
    t0 = 0;
    z0 = [
        nAout0
        nBout0
        nZout0
        Tout0
    ];
    tf = 400; % min
    
    % Solve the ODEs
    [t,zz] = ode45(@odeqns,[t0 tf],z0);
	last_value = length(t);
	t_f = t(last_value);
	z = zz(last_value,:);
	
    % Plot the conversion versus time
    fA = 100*(nA0 - zz(:,1))/nA0;
    figure
    plot(t,fA)
    ylabel('Percent Conversion')
    xlabel('Time (min)')
    
    % also plot T versus time
    figure
    plot(t,zz(:,4))
    ylabel('Temperature (deg R)')
    xlabel('Time (min)')
    	
end % of Activity_23_2.m