% Modified version of the MATLAB template file SolvIVDifD.m used in the
% solution of Example 19.2 of "A First Course on Kinetics and Reaction
% Engineering."
%
function Example_19_2(deltaT)
    % Known quantities and constants (in L, mol, min, J, K, Torr)
    k0 = 265;
    E = 73000;
    dH298 = -165000;
    P0 = 1000;
    T0 = 1225;
    R = 8.314;
    % Constants for heat capacities
    aA = 28;
    aY = 26;
    aZ = 30;
    bA = 0.05;
    bY = 0.01;
    bZ = 0.005;
    V = 1; % basis
    nA0 = P0*V/(62.36367*T0); % Ideal gas law
    Tf = T0 + deltaT;

	% Function that evaluates the ODEs
    function dzdt = odeqns(t,z)
        nA = z(1);
        nY = z(2);
        nZ = z(3);
        T = z(4);
        cpA = aA + bA*T;
        cpY = aY + bY*T;
        cpZ = aZ + bZ*T;
        dH = dH298 + (aZ + 2*aY - 2*aA)*(T-298) + 0.5*(bZ + 2*bY - 2*bA)*(T^2-298^2);
        r = k0*exp(-E/R/T)*(nA/V)^2;
        dzdt = [
           -2*V*r
           2*V*r
           V*r
           -V*r*(dH - R*T)/(nA*(cpA-R) + nY*(cpY-R) + nZ*(cpZ-R));
        ];
    end % of internal function odeqns
        
    % Initial values
    t0 = 0.0;
    z0 = [
        nA0
        0
        0
        T0
    ];

    tf = 100000;
    options = odeset('Events',@stop);
	[t, zz, te, ze, ie] = ode45(@odeqns,[t0, tf],z0,options);
    % NOTE: For stiff equations replace ode45 with ode15s
    
    % Function that provides the integration stopping criterion
    function [stop_when, isterminal, direction] = stop(t,z)
        isterminal = 1;
        direction = 0;
        
        stop_when = Tf - z(4);
    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
    T_f = z(4)
    pct_converted = 100*(nA0 - z(1))/(nA0)
    
end % of Example_19_2.m
