function [nAg,nAl,nBl,PAi] = Example_39_1_bulk(given,guess)
%Example_39_1_bulk Solve the bulk phase mole balances in Example 39.1 of AFCoKaRE

    % Function that evaluates the bulk mole balance equations
    function f = evalEqns(z)
        % Extract the unknowns
        nag = z(1);
        nal = z(2);
        nbl = z(3);
        pai = z(4);
       
        % Calculate quantities appearing in the bulk mole balance equations
        ntotalg = given.P*given.VFRg0/given.R/given.T - given.nAg0 + nag; % eq 10
        VFRg = ntotalg*given.R*given.T/given.P; % eq 11
        PAg = nag/ntotalg*given.P;
        CAl = nal/given.VFRl; % eq 12
        CBl = nbl/given.VFRl; % eq 13
        rate = given.k*given.Vl*CAl*CBl;
        if (rate < 0) 
            rate = 0.;
        end
        
        % Solve the liquid film mole balances
        [NA0, NAxL, NBxL] = Example_39_1_film(given,pai,CAl,CBl);
        
        % Evaluate the bulk mole balance equations
        f = [
            given.nAg0 - given.SV*given.Vl*NA0 - nag; % eq 4
            NAxL*given.SV*given.Vl - rate - nal; % eq 5
            given.nBl0 + NBxL*given.SV*given.Vl - rate - nbl; % eq 6
            given.kGP*(PAg - pai) - NA0; % eq 20
        ];
    end % of internal function evalEqns
    
    % Solve the bulk phase mole balance equations
    z = fsolve(@evalEqns, guess);
    nAg = z(1);
    nAl = z(2);
    nBl = z(3);
    PAi = z(4);
	
end % of Example_39_1_bulk.m