% ME515 Numerical Solns of PDEs, Prof. John Sullivan
% The function accepts as Inputs: A,B,C,R and n
% 	It returns S as the solution;   
% Copy this .m file into the directory containing
% the main program
% Access using the following syntax:
% S = thomas(A,B,C,R,n)
% In the main program.

function S = thomas(A,B,C,R,n)

%  Usage: A,B,C are the three bands of the matrix.  
%	A is left of diagonal
%	B is the diagonal band
%	C is to the right of diagonal
%	R is the "right hand side" vector.
%	The output is S in the form [Matrix]S = {R}
%  Note: The names of the arrays in your main program
%	are not necessarily the same as in the Thomas
%	subroutine.

C(1)=C(1)/B(1);
R(1)=R(1)/B(1);
for i=2:n
    im = i-1;
    C(i)=C(i)/(B(i)-A(i)*C(im));
    R(i)=(R(i)-A(i)*R(im))/(B(i)-A(i)*C(im));
end

for i=1:n-1
    j=n-i;
    jp=j+1;
    R(j)=R(j)-C(j)*R(jp);
end
S=R;