function [H] = hessian(x,t,h,order)
%
%  Hessian of [H3.14]
%
J = jacobian(x,t);
H = eye(size(J));
fac = 1;
JJ = H;
for i=1:order,
  fac = fac*i;
  JJ = JJ*J*h;
  H = H + JJ/fac;
end
return