% MATLAB template file for solving a set of mixed boundary value ordinary
% differential equations of the matrix form
%   dy1/dx = S11*y1/x + ... + S1n*yn/x + f1(x,y1,...,yn)
%   dy2/dx = S21*y1/x + ... + S2n*yn/x + f2(x,y1,...,yn)
%              .
%              .
%              .
%   dyn/dx = Sn1*y1/x + ... + Snn*yn/x + fn(x,y1,...,yn)
% where x ranges from 0 to b and the value of some of the y's are known at
% x = 0 while the values of the other y's are known at x = 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 seven modifications before each use; the locations
% where editing is required are indicated by the comment "% EDIT HERE"
%
function result = SolvBVDifS
    % Known quantities and constants (in consistent units)
% EDIT HERE (Required modification 1 of 7):
    % define universal and problem-specific constants here

    % 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)
% EDIT HERE (Required modification 2 of 7):
		dydx = [
			% Enter expression for f1
			% Enter expression for f2
            % and so on, one fi 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 7):
		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 7):
	% 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 7):
	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
	];

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

% EDIT HERE (Required modification 6 of 7) Provide the matrix, S, in 
% dyi/dx = Si1*y1/x + ... + Sin*yn/x + fn(x,y1,...,yn):
	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);
    
% EDIT HERE (Required modification 7 of 7):
    % 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 SolvBVDifS.m
