function x = hessen(A,b);
% HESSEN  Solve system A x = b by Gauss elimination and
%   back-substitution, without pivoting, on matrices A with
%   upper Hessenberg form    * * * * *
%                            * * * * *
%                       A =  0 * * * *
%                            0 0 * * *
%                            0 0 0 * *

[n m] = size(A);
if n ~= m, error('A not square'), end
if size(b,2) ~= 1, error('b not a column vector'), end
if n ~= size(b), error('b wrong size'), end

for i = 2:n
  for k = i+2:n
    if A(k,i-1) ~= 0, error('not upper Hessenberg form'), end
  end
  r = A(i,i-1) / A(i-1,i-1);             % the one multiplier per row
  A(i,i:n) = A(i,i:n) - r * A(i-1,i:n);  % implied loop over columns
  b(i) = b(i) - r * b(i-1);
end

x = bs1(triu(A),b);            % we've already got a back-sub routine
