function [flux, eta] = Example_38_2_cat_bal(given,CAs)
%Example_38__cat_bal Solve the catalyst phase mole balance for Example 38.2
%of AFCoKaRE and return the rate

    % Structure array given contains known quantities and constants

    % Define the internal functions that will be needed

	% Function that evaluates the column vector f in dydx = S*y/x + f(x,y)
	function dydx = bvodes(x,y)
		dydx = [
            y(2)
            given.k/given.DA*y(1)^2
		];
    end % of internal function bvodes

	% Function that calculates the errors at the boundaries
	function res = bvs(y_at_start,y_at_end)
		res = [
            y_at_start(2)
            y_at_end(1) - CAs
		];

    end % of internal function bvs

    % All the necessary internal functions are now defined

	x_range_low = 0.;
	x_range_high = given.Rp;
	n_mesh_points = 20;
	
	% The next line creates an array for the independent variable
	x = linspace(x_range_low,x_range_high,n_mesh_points);

	% Guesses
	yinit = [
        CAs/2.
        CAs/given.Rp
	];

	% The next line creates a structure containing the mesh and guesses
	solinit=bvpinit(x,yinit);

	S = [
        0,0
        0,-2
	];

    % Solve the ODEs
	options = bvpset('SingularTerm',S);
	soln = bvp4c(@bvodes,@bvs,solinit,options);
    
    lastPoint = length(soln.x);
    flux = -soln.y(2,lastPoint)*given.DA; % solution, equation (8)
    eta = -3*flux/(given.Rp*given.k*CAs^2); % solution, equation (10)

end % of Example_38_2a
