% Modified version of MATLAB template file SolvIVDifI used to model the
% CSTR zone in Activity 35.1
function [nAout,nBout,nZout] = Activity_35_cstr(nAin,nBin,nIin,nZin)
    % Known quantities and constants (in consistent units)
    Ltot = 6.0; % m
    D = 0.07; % m
    V = 0.05*pi()*D^2/4*Ltot; % m^3
    T = 450 + 273.15; % K
    P = 5.0; % atm
    k = 59.5; % mol/h/atm^0.5/m^3
    R = 8.206E-5; % m^3*atm/mol/K
    
    % Function that evaluates the equations
    function f = evalEqns(z)
        nA = z(1);
        nB = z(2);
        nZ = z(3);
        ntot = nA + nB + nZ + nIin;
        CA = nA*P/ntot/R/T;
        CB = nB*P/ntot/R/T;
        r = k*sqrt(CA)*CB;
        f = [
            nA - nAin - 2*V*r
            nB - nBin - V*r
            nZ - nZin + 2*V*r
        ];
    end % of internal function evalEqns
    
    % guesses for the solution
    z_guess = [
        nAin/2
        nBin/2
        nZin + nAin/2
    ];

    % Solve the set of algebraic equations
    z = fsolve(@evalEqns, z_guess);
    display('The solver found the following values for the unknowns:');
    z
    display('The corresponding values of the functions being solved are as follows:');
    f = evalEqns(z)
    
    % Set the return values
    nAout = z(1);
    nBout = z(2);
    nZout = z(3);
	
end % of Activity_35_cstr.m