function ShowHessLU()
clc
format short
% Generate an upper Hessenberg system...
n = 6;
H = triu(randn(n,n),-1)
b = randn(n,1);
% Compute PH = LU...
[m,U,piv] = HessLU(H);
% And use m, U, and piv to solve Hx = b...
x = HessSolve(m,U,piv,b);
%
TheTwoSolutions = [x H\b]

function [m,U,piv] = HessLU(H)
% H is an nxn upper Hessenberg matrix (H(i,j) = 0 if i>j+1)
% This function carries out Gaussian elimination to compute PH = LU
% but encodes P and L through (n-1)-vectors m and piv. In particular,
% M(n-1)P(n-1)***M(1)P(1)H = U (upper triangular) where
% P(i) is the identity if piv(k) = 0. Otherwise, P(k) is the
% identity with rows i and i+1 swapped.
% M(i) is the identity except that its (i+1,i) entry is -m(i).
[n,n] = size(H);
m = zeros(n-1,1);
piv = zeros(n-1,1);
for i=1:n-1
if abs(H(i+1,i))>abs(H(i,i))
piv(i) = 1;
H([i+1 i],i:n) = H([i i+1],i:n);
end
m(i) = H(i+1,i)/H(i,i);
H(i+1,i:n) = H(i+1,i:n) - m(i)*H(i,i:n);
end
U = triu(H);

function x = HessSolve(m,U,piv,b)
% Solves the nxn nonsingular upper Hessenberg system Hx = b where
% m, U, and piv are obtained via [m,U,piv] = HessLU(H);

% First, b <-- M(n-1)P(n-1)****M(1)P(1)b where the swap permutations
% P(i) are encoded in piv and the "single multiplier" ,matrices M(i) are
% encoded in m...
n = length(b);
for i=1:n-1
if piv(i) == 1
b([i+1 i]) = b(i:i+1);
end
b(i+1) = b(i+1)-m(i)*b(i);
end
% Solve Ux = updated b ...
x = U\b;