function x=back_fw_subs(L,U,b)
%   Construct a backward substitution for the linear syste
%      A*x=b
%   where A=L*U and L is lower triangular and U is upper triangular
% 
    
    n=length(b);


    % Compute: z=L\b;
    z=zeros(size(b));
    for k=1:n
        z(k)=(b(k)-L(k,1:(k-1))*z(1:(k-1)))/L(k,k);
    end

    % Compute: x=U\z;
    x=zeros(size(b));
    for k=n:-1:1
        x(k)=  (z(k)- U(k,(k+1):n)*x((k+1):n))/U(k,k);
    end

    