% Modified version of the MATLAB template file SolvIVDifD.m used in the
% solution of Example 19.1 of "A First Course on Kinetics and Reaction
% Engineering."
%
function Example_19_1
    % Known quantities and constants (in L, mol, s, K, J, g)
    k01 = 8110;
    k02 = 9730;
    E1 = 75000;
    E2 = 66000;
    dH1 = -41000;
    dH2 = -36000;
    cp = 1.23*4.184;
    rho = 1020;
    V = 1900;
    T0 = 180 + 273.15;
    CA0 = 2.9;
    CB0 = 3.2;
    CSf = 1.0;
    R = 8.314;

	% Function that evaluates the ODEs
    function dzdt = odeqns(t,z)
        nA = z(1);
        nB = z(2);
        nR = z(3);
        nS = z(4);
        nY = z(5);
        nZ = z(6);
        T = z(7);
        r1 = k01*exp(-E1/R/T)*(nA/V)*(nB/V);
        r2 = k02*exp(-E2/R/T)*(nR/V)*(nB/V);
        dzdt = [
            -V*r1
            V*(-r1 - r2)
            V*(r1 - r2)
            V*r1
            V*r2
            V*r2
            -V*(r1*dH1 + r2*dH2)/(rho*V*cp)
        ];
    end % of internal function odeqns
        
    % Initial values
    t0 = 0;
    z0 = [
        CA0*V
        CB0*V
        0
        0
        0
        0
        T0
    ];

    tf = 100000;
    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;
        

        stop_when = CSf - z(4)/V;
    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,:);
	
    t_f_min = t_f/60
    CS = z(4)/V
    CB = z(2)/V
    CR = z(3)/V
    CZ = z(6)/V
    T_C = z(7)-273.15
    
end % of Example_19_1.m
