% Modified version of the MATLAB template file SolvIVDifI.m used in the
% solution of Example 23.2 of "A First Course on Kinetics and Reaction
% Engineering."
%
function Example_23_2
    % Known quantities and constants (in consistent units)
    V = 500; % cm3
    VFR = 1; % cm3/s
    CAin = 0.015; % mol/cm3
    CBin = CAin;
    Tin = 50 + 273.15; % K
    cp = 0.35*4.1840; % J/g/K
    rho = 0.93; % g/cm3
    dH = -20000; % J/mol
    k0 = 3.24e12; % cm3/mol/s
    E = 25000*4.1840; % J/mol
    R = 8.31446; % J/mol/K
    Tout0 = 190 + 273.15; % K
    Toutf = 260 + 273.15; % K
    nYin = 0; % mol/s
    nZin = 0; % mol/s
    % Other calculated constants
    nAin = CAin*VFR; % mol/s
    nBin = CBin*VFR; % mol/s
    
	% 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);
        % Other variable quantities
        CA = nA/VFR;
        CB = nB/VFR;
        r = k0*exp(-E/(R*T))*CA*CB;
        % Evaluate the ODEs
        dzdt = [
            VFR/V*(nAin - nA - V*r)
            VFR/V*(nBin - nB - V*r)
            VFR/V*(nYin - nY + V*r)
            VFR/V*(nZin - nZ + V*r)
            (-VFR*rho*cp*(T-Tin) - V*r*dH)/(V*rho*cp)
        ];
    end % of internal function odeqns
        
    % Initial values
    t0 = 0;
    z0 = [
        0
        0
        0
        0
        Tout0;
    ];

    tf = 3000; % Don't want to reach this value
    options = odeset('Events',@stop);
	[t, zz, te, ze, ie] = ode45(@odeqns,[t0, tf],z0,options);
    
    % Function that provides the integration stopping criterion
    function [stop_when, isterminal, direction] = stop(t,z)
        isterminal = 1;
        direction = 0;
        % final condition
        stop_when = z(5) - Toutf;
    end % of internal function stop

    % 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,:);
	
    % Report the required time and the conversion
    required_time = t_f
    conversion = (nAin - z(1))/nAin
    
end % of Example_23_2.m
