function x = tridiag(a,b,c,f)
% x = tridiag(a,b,c,f)
%
% Solve an nxn tridiagonal system with sub-diagonal a, diagonal b, 
% superdiagonal c, and rhs f 
%

n = length(f);
x = zeros(size(f));

% Forward elimination
for i=2:n
    p = a(i-1)/b(i-1);
    b(i) = b(i) - p*c(i-1) ;
    f(i) = f(i) - p*f(i-1) ;
end

% Back substitution
x(n) = f(n)/b(n);
for i=n-1:-1:1
   x(i) = (f(i) - c(i)*x(i+1))/b(i);
end