function x = backsub(U, b)
% BACKSUB   Solve an upper triangular linear system.
%
% Input:
%   U    upper triangular square matrix (n by n)
%   b    right-hand side vector (n by 1)   
%
% Output:
%   x    solution of Ux=b (n by 1)
%
    n = length(b);
    x = zeros(n, 1);
    for i = n:-1:1
        x(i) = (b(i) - U(i,i+1:n)*x(i+1:n)) / U(i,i);
    end
end