function reb_9_6_2()
% Calculations for Example 9.6.2 of Reaction Engineering Basics
    % global constants
    V = 10.0E3; % cm^3
    Vex = 1.4E3; % cm^3
    Uex = 138.; % cal /ft^2 /min /K
    Aex = 1200./929.; % ft^2
    Tex0 = 40. + 273.15; % K
    Tex_in = 40. + 273.15; % K
    mDot_ex = 100.; % g /min
    rho = 1.0; % g /cm^3
    rho_ex = 1.0; % g /cm^3
    Cp = 1.0; % cal /g /K
    Cp_ex = 1.0; % cal /g /K
    CA_0 = 5.0E-3; % mol /cm^3
    CB_0 = 7.0E-3; % mol /cm^3
    dH_1 = -16.7E3; % cal /mol
    dH_2 = -14.3E3; % cal /mol
    k0_1 = 9.74E12; % cm^3 /mol /min
    E_1 = 20.1E3; % cal /mol
    k0_2 = 2.38E13; % /min
    E_2 = 25.3E3; % cal /mol
    t_f = 30.; % min
    fA_f = 0.45;
    Re = 1.987; % cal /mol /K

    % residuals function
    function epsilon = residual(guess)
        % extract the individual guesses
        T_0 = guess(1);

        % define the initial values for solving the IVODEs
        ind_0 = 0.0;
        nA_0 = CA_0*V;
        nB_0 = CB_0*V;
        nX_0 = 0.0;
        nY_0 = 0.0;
        nZ_0 = 0.0;
        dep_0 = [nA_0; nB_0; nX_0; nY_0; nZ_0; T_0; Tex0];

        % define the stopping criterion for solving the IVODEs
        stop_var = 0;
        stop_val = t_f;
        
        % solve the IVODEs
        odes_are_stiff = false;
        [ind, dep, flag] = solve_ivodes(ind_0, dep_0, stop_var, stop_val...
            , @derivatives, odes_are_stiff);
    
        % check for solver issues
        if flag <= 0
            disp(' ')
            disp('WARNING: The ODE solution may not be accurate!')
        end

        % extract the calculated final molar amount of A
        nA = dep(:,1);
        nA_f_calc = nA(end);

        % calculate the specified final molar amount of A
        nA_0 = CA_0*V;
        nA_f = nA_0*(1 - fA_f);

        % evaluate the residual
        eps_T0 = nA_f_calc - nA_f;

        % return the residuals as a column vector
        epsilon = [eps_T0];
    end

    % derivatives function
    function ddt = derivatives(ind, dep)
        % extract variables needed to evaluate the derivatives
        nA = dep(1);
        nB = dep(2);
        T = dep(6);
        Tex = dep(7);
        
        % calculate additional unknowns
        % calculate the rates
        CA = nA/V;
        CB = nB/V;
        k_1 = k0_1*exp(-E_1/Re/T);
        r1 = k_1*CA*CB;
        k_2 = k0_2*exp(-E_2/Re/T);
        r2 = k_2*CA;

        % calculate the rate of heat exchange
        Qdot = Uex*Aex*(Tex-T);

        % evaluate the derivatives
        dnAdt = -V*(r1 + r2);
        dnBdt = -V*r1;
        dnXdt = V*r1;
        dnYdt = V*r1;
        dnZdt = V*r2;
        dTdt = (Qdot - V*(r1*dH_1 + r2*dH_2))/rho/V/Cp;
        dTexdt = (-Qdot - mDot_ex*Cp_ex*(Tex-Tex_in))/rho_ex/Vex/Cp_ex;

        % return the derivatives as a column vector
        ddt = [dnAdt; dnBdt; dnXdt; dnYdt; dnZdt; dTdt; dTexdt];
    end

    % BSTR function
    function [T0, t, nA, nB, nX, nY, nZ, T, Tex] = BSTR_variables(T0_guess)
        % define the initial guess for solving the ATEs
        initial_guess = [T0_guess];

        % solve the implicit ATE
        options = optimoptions('fsolve','Display','off');
        [soln, eps, flag, details] = fsolve(@residual, initial_guess, options);
    
        % check for solver issues
        if flag <= 0
            disp(' ')
            disp('WARNING: The ATE solver did not converge')
            disp(['         ' , details.message])
        end

        % extract the solution of the implicit ATE
        T0 = soln(1);

        % define the initial values for solving the IVODEs
        ind_0 = 0.0;
        nA_0 = CA_0*V;
        nB_0 = CB_0*V;
        nX_0 = 0.0;
        nY_0 = 0.0;
        nZ_0 = 0.0;
        dep_0 = [nA_0; nB_0; nX_0; nY_0; nZ_0; T0; Tex0];

        % define the stopping criterion for solving the IVODEs
        stop_var = 0;
        stop_val = t_f;

        % solve the IVODEs
        odes_are_stiff = false;
        [ind, dep, flag] = solve_ivodes(ind_0, dep_0, stop_var, stop_val...
            , @derivatives, odes_are_stiff);
    
        % check for solver issues
        if flag <= 0
            disp(' ')
            disp('WARNING: The ODE solution may not be accurate!')
        end

        % extract the solution of the IVODEs
        t = ind;
        nA = dep(:,1);
        nB = dep(:,2);
        nX = dep(:,3);
        nY = dep(:,4);
        nZ = dep(:,5);
        T = dep(:,6);
        Tex = dep(:,7);
    end

    % deliverables function
    function deliverables()
        % define an initial guess for T0
        T0_guess = Tex_in;

        % solve the BSTR design equations
        [T0, t, nA, nB, nX, nY, nZ, T, Tex] = BSTR_variables(T0_guess);

        % convert T0 to deg C
        T0 = T0 - 273.15;

        % calculate the other quantities of interest
        T_f = T(end) - 273.15;
        Tex = Tex(end) - 273.15;
        sel_X_Z = nX(end)/nZ(end);

        % tabulate the results
        item = ["T0";"T_f";"Te_f";"sel_X_Z"];
        value = [T0; T_f; Tex; sel_X_Z];
        units = ["°C"; "°C"; "°C"; "mol X per mol Z"];
        results_table = table(item,value,units);
    
        % display the results table without the quotes
        str_table = results_table;
        str_table.item = categorical(str_table.item);
        str_table.units = categorical(str_table.units);
        disp(' ')
        disp(str_table)
    
        % save the results
        writetable(results_table,'matlab_results.csv');
    end

    % execution command
    deliverables();
end