% Modified version of the AFCoKaRE MATLAB template file SolvIVDifD.m used
% in the solution of Example 32.2 of "A First Course on Kinetics and
% Reaction Engineering. More specifically, this file performs the
% calculations for the second stage of semibatch operation where the
% pressure relief valve is open.
%
function [ts2,nAs2,nBs2,nYs2,nZs2] = Example_32_2_semi_b(t0,nA0,nB0,nY0,nZ0,tf)
    % Known quantities and constants (in min, mol, L, atm and K)
    K = 8e-5;
    hZ = 8e-4; % mol/L/atm
    k = 5e-4; % L/mol/min
    V = 100*3.785; % L
    CA0 = 4.; % mol/L
    CB0 = 4.; % mol/L
    T = 80 + 273.15; % K
    Vhead = 8*3.785; % L
    Pmax = 1.5; % atm
    % Constants found in handbooks and other reference sources
    R = 0.08206; % L*atm/mol/K

	% Function that evaluates the ODEs
    function dzdt = odeqns(t,z)
        nA = z(1);
        nB = z(2);
        nY = z(3);
        CA = nA/V;
        CB = nB/V;
        CY = nY/V;
        CZ = nZ0/V;
        r = k*CA*CB*(1 - CY*CZ/K/CA/CB);
        dzdt = [
            -V*r
			-V*r
            V*r
        ];
    end % of internal function odeqns
        
    % Initial and final values of the dependent variables (values for the
    % independent variables are passed in as arguments)
    z0 = [
        nA0
        nB0
        nY0
    ];
    
    % Solve the ODEs
    [t,zz] = ode45(@odeqns,[t0 tf],z0);

    % Extract the final values of the independent and dependent variables
	last_value = length(t);
	t_f = t(last_value); % final value of the independent variable
	z = zz(last_value,:); % final values of the dependent variables
	
    % Set the values of the return variables (t is already set)
    ts2 = t;
    nAs2 = zz(:,1);
    nBs2 = zz(:,2);
    nYs2 = zz(:,3);
    nZs2 = nZ0*ones(length(t),1);
    
end % of Example_32_2_semi_b.m
