function x = thomas(a,b,c,rhs)
%thomas Thomas algorithm for solving the tridiagonal system Ax = rhs.

%	 x = thomas(a,b,c,rhs)
%   The function solves a tridiagonal system of linear equations Ax = rhs
%   using the linear Thomas algorithm. a is the lower diagonal, b the
%   diagonal, and c the upper diagonal.

n = length(b);
x = zeros(n,1);

c(1) = c(1)/b(1);
rhs(1) = rhs(1)/b(1);
for i = 2:n-1
    c(i) = c(i)/(b(i)-a(i-1)*c(i-1));
    rhs(i)=(rhs(i)-a(i-1)*rhs(i-1))/(b(i)-a(i-1)*c(i-1));
end
rhs(n)=(rhs(n)-a(n-1)*rhs(n-1))/(b(n)-a(n-1)*c(n-1));

% now perform back substitution
x(n) = rhs(n);
for i = n-1:-1:1
    x(i) = rhs(i)-c(i)*x(i+1);
end
