% Modified version of the MATLAB template file SolvIVDifD.m used in the
% solution of Example 35.2 of "A First Course on Kinetics and Reaction
% Engineering."
%
function Example_35_2
    % Known quantities and constants
    k01 = 8110*60; % L/mol/min
    k02 = 9730*60; % L/mol/min
    E1 = 75000; % J/mol
    E2 = 66000; % J/mol
    V_total = 1900; % L
    T = 190 + 273.15; % K
    CA0 = 2.9; % mol/L
    CB0 = 3.2; % mol/L
    CSf = 1.0; % mol/L
    R = 8.314; % J/mol/K
    VFR = 75.; % L/min
    
    % Calculated constants
    V_R2 = 0.1*V_total; % L
    V_R1 = V_total - V_R2; % L
    k1 = k01*exp(-E1/R/T); % L/mol/min
    k2 = k02*exp(-E2/R/T); % L/mol/min

	% Function that evaluates the ODEs
    function dzdt = odeqns(t,z)
        nA_FS1 = z(1);
        nB_FS1 = z(2);
        nR_FS1 = z(3);
        nS_FS1 = z(4);
        nY_FS1 = z(5);
        nZ_FS1 = z(6);
        nA_FS2 = z(7);
        nB_FS2 = z(8);
        nR_FS2 = z(9);
        nS_FS2 = z(10);
        nY_FS2 = z(11);
        nZ_FS2 = z(12);
        CA_R1 = nA_FS2/VFR;
        CB_R1 = nB_FS2/VFR;
        CR_R1 = nR_FS2/VFR;
        r1_R1 = k1*CA_R1*CB_R1;
        r2_R1 = k2*CR_R1*CB_R1;
        CA_R2 = nA_FS1/VFR;
        CB_R2 = nB_FS1/VFR;
        CR_R2 = nR_FS1/VFR;
        r1_R2 = k1*CA_R2*CB_R2;
        r2_R2 = k2*CR_R2*CB_R2;
        dzdt = [
            VFR/V_R2*(nA_FS2 - nA_FS1 + V_R2*(-r1_R2))
            VFR/V_R2*(nB_FS2 - nB_FS1 + V_R2*(-r1_R2 - r2_R2))
            VFR/V_R2*(nR_FS2 - nR_FS1 + V_R2*(r1_R2 - r2_R2))
            VFR/V_R2*(nS_FS2 - nS_FS1 + V_R2*(r1_R2))
            VFR/V_R2*(nY_FS2 - nY_FS1 + V_R2*(r2_R2))
            VFR/V_R2*(nZ_FS2 - nZ_FS1 + V_R2*(r2_R2))
            VFR/V_R1*(nA_FS1 - nA_FS2 + V_R1*(-r1_R1))
            VFR/V_R1*(nB_FS1 - nB_FS2 + V_R1*(-r1_R1 - r2_R1))
            VFR/V_R1*(nR_FS1 - nR_FS2 + V_R1*(r1_R1 - r2_R1))
            VFR/V_R1*(nS_FS1 - nS_FS2 + V_R1*(r1_R1))
            VFR/V_R1*(nY_FS1 - nY_FS2 + V_R1*(r2_R1))
            VFR/V_R1*(nZ_FS1 - nZ_FS2 + V_R1*(r2_R1))
        ];
    end % of internal function odeqns
        
    % Initial values
    t0 = 0;
    z0 = [
        CA0*VFR
        CB0*VFR
        0
        0
        0
        0
        CA0*VFR
        CB0*VFR
        0
        0
        0
        0
    ];

    tf = 1000;
    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;
        
        nS_FS1 = z(4);
        nS_FS2 = z(10);
        CS_R1 = nS_FS2/VFR;
        CS_R2 = nS_FS1/VFR;
        stop_when = CSf - (V_R1*CS_R1 + V_R2*CS_R2)/V_total;
    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
    
end % of Example_35_2.m
