function hess = powellhess(x)
% usage:  hess = powellhess(x)
% description: Returns the Hessian matrix for
% powell(x).

% a little bit of bulletproofing
if (~isvector(x) | (length(x) ~= 4))
  error('powell: needs a two-dimensional vector argument');
end
if (size(x,2) ~= 1)
  error('powell: needs a column vector');
end 
% the actual calculation
hess = [2+1200*(10*x(1)-x(4))^2, 20, 0, -120*(10*x(1)-x(4))^2;
        20, 200+12*(x(2)-2*x(3))^2, -24*(x(2)-2*x(3))^2, 0;
        0, -24*(x(2)-2*x(3))^2, 10+48*(x(2)-2*x(3))^2, -10; 
        -120*(10*x(1)-x(4))^2, 0, -10, 10+12*(10*x(1)-x(4))^2];