function u=TRIDSOLVE(a,b,c,r,n)
% function that solves for a simultaneous system of equations
% where only the diagonals of P are nonzero
%
% INPUT: a(subdiagonal, 2:n), b(main diagonal), c(supradiagonal, 1:n-1)
% r=right hand side vector, u=vector of unknowns, n=number of unknowns
%
% assumes no zero pivot for the LU i.e. diagonal dominance
%
% Idea from Thomas algorithm, Numerical Recipes
gam=zeros(n,1);
if b(1)==0; fprintf('tridsolve: rewrite equations'); end;
%if this happens, rewrite eqns as set of o(n-1), w/ u2 trivially eliminated
bet=b(1,1);
u(1,1)=r(1,1)/bet;
for j=2:n					% decomposition and forward substitution
   gam(j,1)=c(j-1,1)/bet;
   bet=b(j,1)-a(j,1)*gam(j,1);
   if bet==0; fprintf('tridsolve failed: zero pivot!'); return; end;
   %algorithm failed; zero pivot, and no pivoting - use more gral alg.
   u(j,1)=(r(j,1)-a(j,1)*u(j-1,1))/bet;
end
for j=n-1:-1:1				% backsubstitution
   u(j,1)=u(j,1)-gam(j+1,1)*u(j+1,1);
end
