function x = my_trilu_solve(alpha,beta,c,f)
%
% MY_TRILU_SOLVE - solve tridiag system with MY_TRILU output
%
n = length(c);  x = zeros(n,1);  z = zeros(n,1);

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

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

