function x=trilu_solve(alpha,beta,c,f)
%
%TRILU_SOLVE - solve tridiagonal system using decompostion
%              produced by TRILU.
%
m=length(c);
x=zeros(m,1);
z=zeros(m,1);

% solve Lz=f by forward substitution
z(1)=f(1);
for k=2:m
    z(k)=f(k)-beta(k)*z(k-1);
end

% solve Ux=z by backward substitution
x(m)=z(m)/alpha(m);
for k=m-1:-1:1
    x(k)=(z(k)-c(k)*x(k+1))/alpha(k);
end

