function [NA0, NAxL, NBxL] = Example_39_1_film(given,PAi,CAl,CBl)
%Example_39_1_film Solve the liquid film mole balances in Example 39.1 of AFCoKaRE

	% Function that evaluates the derivatives
	function dydx = bvodes(x,y)
        rate = given.k*y(1)*y(3);
        if (rate < 0) 
            rate = 0;
        end
		dydx = [
            y(2)
            rate/given.DA
            y(4)
            rate/given.DB
		];
    end % of internal function bvodes

	% Function that calculates the errors at the boundaries
	function res = bvs(y_at_start,y_at_end)
        % y_at_start is a column vector containing y values at the starting
        % x boundary and y_at_end is a column vector containing y values at
        % the ending x boundary.
		res = [
			y_at_start(1) - PAi/given.hA
            y_at_end(1) - CAl
            y_at_start(4)
            y_at_end(3) - CBl
		];

    end % of internal function bvs

	% Set up the initial mesh
	x_range_low = 0.0;
	x_range_high = given.xL;
	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 = [
		(PAi/given.hA + CAl)/2
        (PAi/given.hA - CAl)/given.xL
        CBl/2
        -CBl/given.xL
	];

	% Create a structure containing the mesh and guesses
	solinit=bvpinit(x,yinit);

	% Solve the odes
	result = bvp4c(@bvodes,@bvs,solinit);
    
    % Calculate any other desired quantities from the results, noting
	% result.x(i) contains the final x value of mesh point i
	% result.y(i,j) contains the final value of yj at mesh point i
    i_last_mesh = length(result.x);
    NA0 = -given.DA*result.y(2,1);
    NAxL = -given.DA*result.y(2,i_last_mesh);
    NBxL = -given.DB*result.y(4,i_last_mesh);

end % of file Example_39_1_film.m
