echo on
%%%%% tridiag_fact.m %%%%%
% Run script and study the effects of each line.
%-From standard Gaussian elimination derive a routine for solving tridiagonal 
% system. 
%-Then include partial pivoting. 
%-Adapt the program also for symmetric positive definite systems using 
% Cholesky factorization.

format compact

n=5;	%1-D system dimension
bc=[1 0 0 0 1]   % boundary condition - right-hand side
% Tridiagonal matrix saved by diagonals
b=[4 4 4 4 4]'; %main diagonal
a=[0 -1 -1 -1 -1]';
c=[-1 -1 -1 -1]';
A=diag(a(2:5),-1)+diag(b,0)+diag(c,1)
x=A\b	%The MATLAB solution

d=zeros(n,1);
m=zeros(n,1);

%%%% Gaussian elimination for tridiagonal system
echo off;
d(1)=b(1);	% diagonal of U
for k=2:n		% for all columns (only one remaining nonzero row)
      m(k)=a(k)/d(k-1);	% calculate multiplier (subdiagonal of L)
      d(k)=b(k) - m(k)*c(k-1); % apply transformation
end
echo on;
% Resulting triangular factors
m	% subdiagonal of L 
d	% diagonal of U
c	% subdiagonal of U is the same as of the original matrix

% If Pivoting is necessary, then the largest element in the column (b or a) has to be chosen
% for pivot.

% To get the solution x, proceed by forward and backward substitution of LUx=b
% for two systems: Ly=b and Ux=y.
% Construct L and U
L=eye(5)+diag(m(2:5),-1)
U=diag(d,0)+diag(c,1)
y=L\b
x=U\y		%Solution again