% Modified version of the MATLAB template file SolvBVDif.m that was used to
% solve Example 1 of Supplemental Unit S6 of "A First Course on Kinetics
% and Reaction Engineering."
%
function result = S6_Example_1
    % Known quantities and constants (in consistent units)
	D = 8.0E-6;
	us = 0.01;
	k = 0.012;
	K = 1.0;
	CA0 = 1.0;
	L = 1.25;

    % Define the internal functions that will be needed

	% Function that evaluates the derivatives
	function dydx = bvodes(x,y)
        dydx = [
            y(2)
            (1/D)*((k+k/K)*y(1)+us*y(2)-k*CA0/K)
		];
    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 = [
            us*y_at_start(1)-D*y_at_start(2)-us*CA0
            y_at_end(2)
        ];

    end % of internal function bvs

    % All the necessary internal functions are now defined

	% Set up the initial mesh
	x_range_low = 0.0;
	x_range_high = L;
	n_mesh_points = 40;
	
	% The next line creates an array for the independent variable
	x = linspace(x_range_low,x_range_high,n_mesh_points);

	% Guesses
	yinit = [
		CA0
		-CA0/L
	];

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

	% Solve the odes
	result = bvp4c(@bvodes,@bvs,solinit);
    
    % Calculate the value of CA (y(1)) at z (x) = L
	lastPoint = length(result.x);
	CAf = result.y(1,lastPoint)
	% Plot the concentration as a function of axial position
	plot(result.x,result.y(1,:))

end % of file S6_Example_1.m
