% Modified version of the MATLAB template file SolvNonDif.m used in the
% solution of Activity 22.2 of "A First Course on Kinetics and Reaction
% Engineering."
%
function [z] = Activity_22_2a(tau)
    % Known quantities and constants (cal, mol, K, min, ft2, kg, L)
    k0 = 2.59e9; % min-1
    E = 16500; % cal/mol
    R = 1.987; % cal/mol/K
    dH = -22200; % cal/mol
    V = 4; % L
    CA0 = 2; % mol/l
    Te0 = 20 + 273.15; % K
    cp = 440; % cal/L/K
    Ve = 0.5; %L
    A = 0.6; % ft2
    U = 1.13e4/60; % cal/ft2/min/K
    rho = 1; % kg/L
    cpe = 1000; % cal/kg/K
    m = 0.2; % kg/min
    T0 = 23 + 273.15; % K
    % Calculated constants
    VFR = V/tau; % L/min
    nA0 = CA0*VFR; % mol/min
    nB0 = 0; % mol/min
    
    % Function that evaluates the equations
    function f = evalEqns(z)
        % mapping of unknowns into vector z
        nA = z(1);
        nB = z(2);
        T = z(3);
        Te = z(4);
        % calculate variables quantities
        CA = nA/VFR;
        Q = U*A*(Te-T);
        r1 = k0*exp(-E/R/T)*CA;
        f = [
            nA0 - nA - V*r1
            nB0 - nB + V*r1
            VFR*cp*(T - T0) + V*r1*dH - Q
            m*cpe*(Te0-Te)-Q
        ];
    end % of internal function evalEqns
    
    % guesses for the solution
    z_guess = [
        nA0
        0
        T0
        Te0
    ];

    % Solve the set of algebraic equations and display the results
    z = fsolve(@evalEqns, z_guess);
    fA = (nA0 - z(1))/nA0;
    nB = z(2);
    	
end % of Example_22_2a.m