function x = backsub(U,z)
% BACKSUB   Forward substitution for lower-triangular linear systems.
% Input:
%   U    upper triangular square matrix (n by n)
%   z    right-hand side vector (n by 1)   
% Output:
%   x    solution of Ux=z (n by 1 vector)

n = length(U);
x = zeros(n,1);
for i = n:-1:1
  x(i) = ( z(i) - U(i,i:n)*x(i:n) ) / U(i,i);
end