% Function to solve the design equations for Example 20.2 from "A First
% Course on Kinetics and Reaction Engineering" for both steps of the
% operating protocol
function Example_20_2
    % Define variables for plotting
    m = linspace(0.1,2.0);
    r_net = zeros(100,1);
    mmax = 0;
    rmax = 0;
    Tmax = 0;
    
    % Calculate the net rate for each coolant flow rate
    for i=1:100
        [tf1,zf1] = Example_20_2_a;
        [tf2,zf2,Tm] = Example_20_2_b(m(i),zf1);
        r_net(i) = zf2(2)/4/(tf1 + tf2 + 25);
        if(r_net(i) > rmax)
            rmax = r_net(i);
            mmax = m(i);
            Tmax = Tm;
        end
    end
    
    % Report optimum flow rate
    display(['Optimum coolant flow (kg/min): ',num2str(mmax)])

    % Report maximum rate
    display(['Maximum net rate of B production (mol/min): ',num2str(rmax)])
    
    % Report maximum reaction temperature
    display(['Maximum reaction temperature (deg C) at optimum flow: ',num2str(Tmax-273.15)])
    
    % Plot the results
    figure
    plot(m,r_net)
    xlabel('Coolant Flow Rate (kg/min)')
    ylabel('Net Rate of B production (mol/min)')

end % of Example_20_2

