function u = tridiag(a, b, c, r, n)
% TRIDIAG:  u = tridiag(a, b, c, r, n)
%   Solves linear tridiagonal systems by Gauss elimination
%   without pivoting.  An example code--A\b surely does it better. 
% (Ed Bueler; 2/27/02)
%
% Returns u which solves A u = r, where
%   b is the diag, c the super diag, and a the sub diag   of A. 
%   Here A is n X n, so a, b, c, u, r are all n x 1 row vects.
%   Note a(1) and c(n) are irrelevant.  They do not 
%   correspond to matrix elements of A.

if b(1)==0.0, error('Error in tridiag: b(1) = 0'), end;
bet=b(1);
u(1)=r(1)/bet;
for j = 2:n
   gam(j)=c(j-1)/bet;
   bet=b(j)-a(j)*gam(j);
   if bet==0.0, error('Error in tridiag: zero in diagonal'), end;
   u(j)=(r(j)-a(j)*u(j-1))/bet;
end;
for j = n-1:-1:1
   u(j)=u(j)-gam(j+1)*u(j+1);
end;