% Modified version of the MATLAB template file SolvBVDifS.m used to solve
% Example 2 of Supplemental Unit S6 of "A First Course on Kinetics and
% Reaction Engineering."
%
function result = S6_Example_2
    % Known quantities and constants (in consistent units)
	Dea = 3.0E-7;
	k = 0.02;
	Rp = 0.3;
	CAs = 1.0;

    % 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)
        	k/Dea*y(1)
		];
    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(2)
			y_at_end(1) - CAs
		];
    end % of internal function bvs

    % All the necessary internal functions are now defined

	% Set up the initial mesh
	x_range_low = 0;
	x_range_high = 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
		CAs/Rp
	];

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

	S = [
	    % Enter comma-separated values S(1,1), S(1,2), ... S(1,N)
	    % Enter S(2,1), S(2,2), ... S(2,N)
	    % Continue entering rows up to S(N,1), S(N,2), ... S(N,N)
	];

    % Solve the ODEs
	options = bvpset('SingularTerm',S);
	result = bvp4c(@bvodes,@bvs,solinit,options);
    
	% Calculate the steady state rate:
	lastPoint = length(result.x);
	rate = (-3*Dea/Rp)*result.y(2,lastPoint)

end % of file S6_Example_2.m
