function x = trisolve(L, U, c, rhs)
% trisolve Solve the equation Tx = b, where T is a tridiagonal matrix.
%
%   x = trisolve(L,U,c,rhs) solves the tridiagonal system Tx = b where T
%   has been factored into a unit lower bidiagonal matrix times an upper
%   bidiagonal matrix. L is the subdiagonal of the lower diagonal matrix,
%   and U is the diagonal of the upper diagonal matrix. c is the
%   superdiagonal of the original tridiagonal matrix.

n = length(U);
x = zeros(n,1);
y = zeros(n,1);

%forward substitution
y(1) = rhs(1);
for i=2:n
	y(i) = rhs(i) - L(i-1)*y(i-1);   
end

%back substitution
x(n) = y(n)/U(n);
for i = n-1:-1:1       
	x(i) =(y(i) - c(i)*x(i+1))/U(i); 
end
end