{ ******************************************************************
  Numerical hessian and gradient
  ****************************************************************** }

procedure HessGrad(X, G : TVector; H : TMatrix);

const
  Eta = 1.0E-6;     { Relative increment }

var
  Delta, Xminus, Xplus, Fminus, Fplus : TVector;
  Temp1, Temp2, F, F2plus             : Float;
  I, J                                : Integer;

begin
  DimVector(Delta, Nvar);   { Increments   }
  DimVector(Xminus, Nvar);  { X - Delta    }
  DimVector(Xplus, Nvar);   { X + Delta    }
  DimVector(Fminus, Nvar);  { F(X - Delta) }
  DimVector(Fplus, Nvar);   { F(X + Delta) }

  F := Func(X);

  for I := 1 to Nvar do
    begin
      if X[I] <> 0.0 then
        Delta[I] := Eta * Abs(X[I])
      else
        Delta[I] := Eta;
      Xplus[I] := X[I] + Delta[I];
      Xminus[I] := X[I] - Delta[I];
    end;

  for I := 1 to Nvar do
    begin
      Temp1 := X[I];
      X[I] := Xminus[I];
      Fminus[I] := Func(X);
      X[I] := Xplus[I];
      Fplus[I] := Func(X);
      X[I] := Temp1;
    end;

  for I := 1 to Nvar do
    begin
      G[I] := (Fplus[I] - Fminus[I]) / (2.0 * Delta[I]);
      H[I,I] := (Fplus[I] + Fminus[I] - 2.0 * F) / Sqr(Delta[I]);
    end;

  for I := 1 to Pred(Nvar) do
    begin
      Temp1 := X[I];
      X[I] := Xplus[I];
      for J := Succ(I) to Nvar do
        begin
          Temp2 := X[J];
          X[J] := Xplus[J];
          F2plus := Func(X);
          H[I,J] := (F2plus - Fplus[I] - Fplus[J] + F) / (Delta[I] * Delta[J]);
          H[J,I] := H[I,J];
          X[J] := Temp2;
        end;
      X[I] := Temp1;
    end;

end;