function U=tridag(a,b,c,r);
% 
%    | b1  c1  0	...	|
%    | a2  b2  c2  0   		|
% A= | 0   a3  b3  c3 0		|;	A*U=r; 
%    | ...			|	Solve for U when A and r are given,
%    | 0   0		 aN  bN |	and A is tridiagonal.
if nargin<4
	disp('needs a, b, c and r as inputs.')
end
n=length(b);
if length(a)==n-1
	a=[0 a];
elseif length(a)==n
	a(1)=0;
end
if length(b)~=length(r)
	disp('forcing vector must have same dimesion as matrix.')
end
if length(c)<n-1
	disp('c vector too short.')
end
% Numerical Recipes, Press, Flannery, Teukolsky and Vetterling,
% Cambridge Univerisity Press (1986), section 2.6, page 40.
% Simple LU decomposition algorithm. Very pretty!
bet=b(1);
if bet==0 
	disp('b1 cannot be zero ...');
	bet=eps;
end
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
	disp('... probably need pivoting.');
	bet=eps;
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;