% tridiagSolve(l,d,u,b)
%
% solves the matrix equation Ax=b via Gaussian elimination...
% where A is the tridiagonal matrix with d[] on the diagonal,
% l[] on the lower-diagonal, and u[] on the upper-diagonal
% d has n-elements, and l and u each have n-1
function x=tridiagSolve(l, d, u, b)
    n=length(d);
    
    % zero below the diagonal
    for i=2:n
        ratio=l(i-1)/d(i-1);
        d(i)=d(i)-ratio*u(i-1);
        b(i)=b(i)-ratio*b(i-1);
    end
    
    x=zeros(n,1);
    x(n)=b(n)/d(n);
    
    % back-substitute to solve
    for i=n-1:-1:1
        x(i)=(b(i)-u(i)*x(i+1))/d(i);
    end
    x = x;
end