function [x]=backsub(U,rhs)
% performs back substitution 
% given an upper triangular matrix U and a rhs vector
% outputs the solution of U*x=rhs

N=length(rhs);

% initialize x as a column vector of the correct length
x=zeros(N,1);

% last row gives last component of x
x(N)=rhs(N)/U(N,N);

for k=N-1:-1:1
    x(k)=rhs(k);
    for j=k+1:N % or j=N:-1:k+1
    % repeated subtraction of coefficients times x-components already computed
    x(k)=x(k)-U(k,j)*x(j);
    end
    x(k)=x(k)/U(k,k);
end

end
    