  function x = HessSolve(alpha,piv,U,b)
% Solves Ax = b where A is an upper Hessenberg whose LU factorization
% PA = LU is encoded (via HessLU) in alpha,  piv, and U.

n = length(b);
for k=1:n-1
    %interchange equations k and k+1 if necessary...
    if piv(k)==1
        b([k+1 k]) = b([k k+1]);
    end
    % apply kth multiplier matrix to b...
    b(k+1) = b(k+1) - alpha(k)*b(k);
end
% Solve Ux = transformed b ... 
x = U\b;