function x = tridiag(e,f,g,b)
% Tridiag: Tridiagonal equation solver for banded system
% Slightly modified from Chapra's code by LTR: see input notes below
% x = tridiag(e,f,g,r): Tridiagonal system solver.
% input:
% e = subdiagonal vector (length n-1)
% f = diagonal vector (length n)
% g = superdiagonal vector (length n-1)
% b = right hand side vector (length n)
% output:
% x = solution vector (column vector, length n)
n=length(f);
if n ~= length(b) || n ~= length(e)+1 || length(e) ~= length(g)
    error('Input vectors are not the correct length.')
end

% Pad the inputs e and g so they're the right length for the tridiag
% algorithm.
try
    e = [0 e];
catch
    e = [0 e'];
end
g(end+1) = 0;

% forward elimination
for k = 2:n
    factor = e(k)/f(k-1);
    f(k) = f(k) - factor*g(k-1);
    b(k) = b(k) - factor*b(k-1);
end
% back substitution
x = zeros(n,1);
x(n) = b(n)/f(n);
for k = n-1:-1:1
    x(k) = (b(k)-g(k)*x(k+1))/f(k);
end