function x = tridiag(DL, D, DU,b)
%  This routine determines the solution of a tridiagonal system of equations.
%  The matrix is listed in three column vectors:
%       L = lower diagonal entries...ie. terms of the form a(i,i-1)
%       D = diagonal entries...ie. terms of the form a(i,i)
%       U = upper diagonal entries...ie. terms of the form a(i,i+1)
%  The right-hand-side is given by the vector "b" and the solution is stored in the vector "x"

n = max(size(b));

L(1) = D(1);

U(1) = DU(1)/L(1);

for k=2:n-1
  L(k) = D(k) - DL(k)*U(k-1);
  U(k) = DU(k)/L(k);
end

L(n) = D(n) - DL(n)*U(n-1);

x(1) = b(1)/L(1);
for k=2:n
  x(k) = (b(k) - DL(k)*x(k-1))/L(k);
end
for k=n-1:-1:1
  x(k) = x(k) - U(k)*x(k+1);
end
