function [ x ] = Usolve( U,y )
%Usolve: solve Ux = y (for x) where U is upper triangular
% note: it is not checked that U is actually upper triangular
    [yRows,yCols] = size(y);
    [URows,UCols] = size(U);
 
    if URows ~= UCols || yRows ~= URows || yCols ~= 1
        error('Either L or b is the wrong size.')
    end
    if any( diag(U) ==  0 )
        error('There are zeros on the diagonal of U')
    end
    
    n = URows;
    x = y;
    for k = n:-1:1  % loop downwards from n to 1
      for j = k+1:n
        x(k) = x(k) - U(k,j)*x(j);
      end
      x(k) = x(k)/U(k,k);
    end
end

