% Modified version of the AFCoKaRE MATLAB template file SolvNonDif.m used
% in the solution of Activity 28.1 of "A First Course in Kinetics and
% Reaction Engineering."
%
function [conv,vol,rate] = Activity_28_1_cstr
    % Known quantities and constants
    k0 = 1E8; % /min
    E = 50000; % J/mol
    dH = -180000; % J/mol
    cP = 3950; % J/l/K
    T0 = 298.; % K
    CA0 = 1.0; % mol/L
    conv = linspace(.2,.9,100)';
    nZ0 = 0; % mol/min
    VFR = 300000/24/60; % L/min
    % Universal constants
    R = 8.3145; % J/mol/K
    % Calculated constant
    nA0 = CA0*VFR;
    
    % Storage for volumes and rates
    vol = zeros(100,1);
    rate = zeros(100,1);
    
    V_last = 100;
    nZ_last = 0.2*nA0;
    T_last = T0 + 23.;
    % Loop through the desired conversions
    fA = 0.2;
    for i = 1:100
        fA = conv(i);
        % guesses for the solution
        z_guess = [
            V_last
            nZ_last
            T_last
        ];

        % Solve the set of algebraic equations
        options = optimoptions(@fsolve,'MaxFunEvals',1000,'Display','off');
        z = fsolve(@evalEqns, z_guess,options);

        % Save the results
        vol(i) = z(1);
        V_last = z(1);
        nZ_last = z(2);
        T_last = z(3);
        rate(i) = k0*exp(-E/R/T_last)*nA0*(1-fA)/VFR;
    end
    
    % Function that evaluates the equations
    function f = evalEqns(z)
        V = z(1);
        nZ = z(2);
        T = z(3);
        nA = nA0*(1-fA);
        CA = nA/VFR;
        r = k0*exp(-E/R/T)*CA;
        f = [
            nA0 - nA - V*r
            nZ0 - nZ + V*r
            VFR*cP*(T-T0) + V*r*dH
        ];
    end % of internal function evalEqns
	
end % of Activity_28_1_cstr