Procedure DFPMIN(var P:array of real; N:integer; FTOL:real;
                                  var ITER:integer; var FRET:real);
const
    ITMAX = 200;  EPS = 0.1e-9;
var
    HESSIN:matrx2;  XI,G,DG,HDG:array[0..50] of real;
    I,J,ITS:integer;    FP,FAC,FAE,FAD,AAA,BBB:real;
begin
    SetLength(HESSIN,51,51);
    FP:=FUNC2(P, N);
    DFUNC(P, G);
    For I:=1 To N do
    begin
        For J:=1 To N do
            HESSIN[I, J]:=0; 
        HESSIN[I, I]:=1; 
        XI[I]:=-G[I];
    end; 
    For ITS:=1 To ITMAX do
    begin
        ITER:=ITS;
        LINMIN(P, XI, N, FRET);
        If 2 * Abs(FRET - FP) <= FTOL * (Abs(FRET) + Abs(FP) + EPS) Then
              Exit;
        FP:=FRET;
        For I:=1 To N do
            DG[I]:=G[I];
        FRET:=FUNC2(P, N);
        DFUNC(P, G);
        For I:=1 To N do
            DG[I]:=G[I] - DG[I];
        For I:=1 To N do
        begin
            HDG[I]:=0;
            For J:=1 To N do
                HDG[I]:=HDG[I] + HESSIN[I, J] * DG[J];
        end; 
        FAC:=0; 
        FAE:=0;
        For I:=1 To N do
        begin
            FAC:=FAC + DG[I] * XI[I];
            FAE:=FAE + DG[I] * HDG[I];
        end; 
        FAC:=1 / FAC;
        FAD:=1 / FAE;
        For I:=1 To N do
            DG[I]:=FAC * XI[I] - FAD * HDG[I];
        For I:=1 To N do
        begin
            For J:=1 To N do
            begin
              AAA:=FAC * XI[I] * XI[J] - FAD * HDG[I] * HDG[J];
              BBB:=FAE * DG[I] * DG[J];
              HESSIN[I, J]:=HESSIN[I, J] + AAA + BBB;
            end; 
        end; 
        For I:=1 To N do
        begin
            XI[I]:=0; 
            For J:=1 To N do
                XI[I]:=XI[I] - HESSIN[I, J] * G[J];
        end; 
    end; 
    If ITS > ITMAX Then ShowMessage('too many iterations in DFPMIN');
end;