function x=supsolve(U,b);
n=size(U,1);
x=zeros(n,1);
for k=n:-1:1
    y=b(k);
    for j=k+1:n
        y=y-U(k,j)*x(j);
    end
    x(k)=y/U(k,k);
end