% Modified version of the MATLAB template file SolvIVDifD.m used in the
% solution of Example 19.3 of "A First Course on Kinetics and Reaction
% Engineering."
%
function Example_19_3(m)
    % 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
    CA0_step1 = 2; % mol/l
    nA0 = 0.82*CA0_step1*V; % mol
    nB0 = 0.18*CA0_step1*V; % mol
    Te0 = 20 + 273.15; %K
    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
    T0 = 60 + 273.15; %K
    CAf = 0.2; % mol/L
    % m in kg/min is passed in as an argument
    

	% 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;
    z0 = [
        nA0
        nB0
        T0
        T0
    ];
    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 = CAf - z(1)/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,:);
	
    % Plot CA vs. t
    CA = zz(:,1)/V;
    figure
    plot(t,CA)
    title(['Coolant flow = ',num2str(m),' kg/min'])
    xlabel('t (min)')
    ylabel('CA (M)')
    % Plot Te and T vs. t
    figure
    plot(t,zz(:,3),t,zz(:,4))
    title(['Coolant flow = ',num2str(m),' kg/min'])
    xlabel('t (min)')
    ylabel('T (K)')
    legend('Reactor','Jacket')
    
end % of Example_19_3
