% MATLAB template file for solving a set of boundary value first order ordinary
% ordinary differential equations of the form
%   dy1/dx = f1(x,y1,...,yn)
%   dy2/dx = f2(x,y1,...,yn)
%              .
%              .
%              .
%   dyn/dx = fn(x,y1,...,yn)
% where x ranges from a to b and the value of some of the y's are known at
% a while the values of the other y's are known at b. 
%
% See Supplemental Unit S9 of "A First Course on Kinetics and Reaction
% Engineering" and/or the MATLAB documentation for additional information
% or if questions arise.
%
% Carl R. F. Lund, 4/22/13
%   last rev. 6/19/13
%
% Returned quantities
%   ans (structure) with the following members
%       ans.x(i) (vector) - value of x at mesh point i
%       ans.y(i,j) (matrix) - value of yj at mesh point i
%
% This file requires six modifications before each use; the locations
% where editing is required are indicated by the comment "% EDIT HERE"
%
function result = SolvBVDif
    % Known quantities and constants (in consistent units)
% EDIT HERE (Required modification 1 of 6):
    % define universal and problem-specific constants here

    % Define the internal functions that will be needed

	% Function that evaluates the derivatives
	function dydx = bvodes(x,y)
% EDIT HERE (Required modification 2 of 6):
		dydx = [
			% Evaluate dy1/dx here
			% Evaluate dy2/dx here
            % and so on, one dyi/dx per line
		];
    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.
% EDIT HERE (Required modification 3 of 6):
		res = [
			% Calculate the error in the first boundary condition here
			% Calculate the error in second boundary condition here
			% and so on, one boundary condition per line
		];

    end % of internal function bvs

    % All the necessary internal functions are now defined

% EDIT HERE (Required modification 4 of 6):
	% Set up the initial mesh
	x_range_low = % insert starting value of the valid range of x here
	x_range_high = % insert ending value of the valid range of x here
	n_mesh_points = % insert the number of mesh points here
	
	% The next line creates an array for the independent variable
	x = linspace(x_range_low,x_range_high,n_mesh_points);

	% Guesses
% EDIT HERE (Required modification 5 of 6):
	yinit = [
		% Enter a guess for y1 (one value for the whole x range) here
		% Enter a guess for y2 (one value for the whole x range) here
        % and so on, one guess per line
	];

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

	% Solve the odes
	result = bvp4c(@bvodes,@bvs,solinit);
    
% EDIT HERE (Required modification 6 of 6):
    % 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

end % of file SolvBVDif.m
