% Modified version of the MATLAB template file SolvIVDifD.m used in the
% solution of Example 20.1 of "A First Course on Kinetics and Reaction
% Engineering." This function models the second step of the operational
% protocol described in that example.
%
function [t_f,z] = Example_20_1_b(z0)
    % Argument z0 contains the initial values of the dependent variables
    % from execution of Example_20_1_a
    
    % Known quantities and constants (cal, mol, K, min, ft2, kg, L)
    k0 = 2.59e9; % min-1
    E = 16500; % cal/mol
    R = 1.987; % cal/mol/K
    dH = -22200; % cal/mol
    V = 4; % L
    cp = 440; % cal/L/K
    Ve = 0.5; %L
    A = 0.6; % ft2
    U = 1.13e4/60; % cal/ft2/min/K
    rho = 1; % kg/L
    cpe = 1000; % cal/kg/K
    Tf = 25 + 273.15; % K
    m = 0.5; % kg/min
    Te0 = 20. + 273.15;

	% Function that evaluates the ODEs
    function dzdt = odeqns(t,z)
        nA = z(1);
        nB = z(2);
        T = z(3);
        Te = z(4);
        r1 = k0*exp(-E/R/T)*nA/V;
        dzdt = [
            -r1*V
            r1*V
            (U*A*(Te-T)-V*dH*r1)/V/cp
            (m*cpe*(Te0-Te)-U*A*(Te-T))/rho/Ve/cpe
        ];
    end % of internal function odeqns
        
    % Initial values
    t0 = 0.0;
    tf = 5000.0;
    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 = Tf - z(3);
    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,:);
	
    % Plot CA vs. t
    CA = zz(:,1)/V;
    figure
    plot(t,CA)
    title(['Cooling Step'])
    xlabel('t (min)')
    ylabel('CA (M)')
    % Plot Te and T vs. t
    figure
    plot(t,zz(:,3),t,zz(:,4))
    title(['Cooling Step'])
    xlabel('t (min)')
    ylabel('T (K)')
    legend('Reactor','Jacket')
    
end % of Example_20_1_b
