function [NA0, NAxL, NBxL] = Activity_39_1_film(given,CAl)
%Activity_39_1_film Solve the liquid film mole balances in Activity 39.1 of AFCoKaRE

	% Function that evaluates the derivatives
	function dydx = bvodes(x,y)
		dydx = [
            y(2)
            given.k/given.DA*y(1)*y(3)
            y(4)
            given.k/given.DB*y(1)*y(3)
		];
    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(1) - given.CAi
            y_at_end(1) - CAl
            y_at_start(4)
            y_at_end(3) - given.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 = [
		(given.CAi + CAl)/2
        (given.CAi - CAl)/given.xL
        given.CBl/10
        -(given.CAi - CAl)/given.xL
	];

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

	% Solve the odes
	result = bvp4c(@bvodes,@bvs,solinit);
    
    % Calculate the return values
    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 Activity_39_1_film.m
