% Modified version of the MATLAB template file SolvIVDifI.m used in the
% solution of Example 23.1 of "A First Course on Kinetics and Reaction
% Engineering."
%
function Example_23_1(VFR, mXout0, mSout0)
    % Known quantities and constants (in consistent units)
    CS0 = 0.04; % g/cm3
    mSin = VFR*CS0; % g/min
    mXin = 0; % g/min
    V = 4.430e3; % cm3
    
	% Function that evaluates the ODEs
    function dzdt = odeqns(t,z)
        mXout = z(1);
        mSout = z(2);
        CX = mXout/VFR;
        CS = mSout/VFR;
        r = 0.014*CS*CX/(0.001 + CS);
        dzdt = [
            VFR/V*(mXin - mXout + V*r)
            VFR/V*(mSin - mSout - 2.2*V*r)
        ];
    end % of internal function odeqns
        
    % Initial and final values
    t0 = 0;
    z0 = [
        mXout0
        mSout0
    ];
    tf = 50000; % min
    
    % Solve the ODEs
    [t,zz] = ode45(@odeqns,[t0 tf],z0);

    % 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 the outlet cell mass flow versus time
    figure
    plot(t,zz(:,1))
    title(['Inlet Flow: ',num2str(VFR),' cm3/min'])
    xlabel('time (min)')
    ylabel('Outlet Cell Mass Flow (g/min)')
    
end % of Example_23_1.m