Метод эверхарта Опции

[-Александр-]

Подскажите где можно почитать по подробней, и посмотреть чью нибудь реализацию все равно ни каком языке.
А то попытался сам реализовать, вроде все сделал как надо а работает хуже метода Эйлера 1 первого порядка.
Вот часть кода, если так с первого взгляда не найдете ошибок, то выложу весь код:
П.С. переменная Step количество учитываемых коэфицентов в полиноме Эверхарта
В комментариях находиться алгоритм выбора шага, он не работает.

Procedure DIFFUN1;
begin
for i:=0 to Nn-1 do for j:=1 to dim do D[i*dim+j]:=0.0;
for i:=0 to Nn-1 do
for ii:=0 to Nn-1 do
if i<>ii then
if C=2 then begin
if I+1<N0 then F1 else
if II+1<N0 then F1; end
else F1;
for i:=0 to Nn-1 do
for j:=1 to dim do
D[i*dim+J+N05]:=Yy[i*dim+j];
end;

Function SQRR(Cg:extended; St:longint):extended;
var Pi:longint; Stt:extended;
begin
Stt:=Cg;
for Pi:=1 to St-1 do Cg:=Cg*Stt;
if St=0 then Cg:=1;
if (St<0) then Cg:=1/Cg;
SQRR:=Cg;
end;

Function XYHK(Hj:extended):extended;
var pi:longint;
begin
for i:=1 to N do begin
YY[i]:=Y0[i]+D[I]*Hj;
for Pi:=1 to Step-1 do
YY[i]:=YY[i]+(A[i,Pi]*SQRR(Hj,Pi+1))/(Pi+1);
end;
end;

Function Ava(Lk:longint):extended;
var Pi,Pj:longint;
begin
for i:=1 to N do
for Pi:=1 to Lk do begin
A[i,Pi]:=Aa[i,Pi];
for Pj:=Pi to Step-2 do
A[i,Pi]:=A[i,Pi]+Aa[i,Pj+1]*Ca[Pj+1,Pi]; end;
end;

Function Ama(Lk:longint):extended;
var Pi,Pj:longint;
begin
for i:=1 to N do begin
Aa[i,Lk]:=(D[i]-D0[i])/Thk[Lk+1];
for Pi:=1 to Lk-1 do
Aa[i,Lk]:=(Aa[i,Lk]-Aa[i,Pi])/(Thk[Lk+1]-Thk[Pi+1]);
end;
end;

Procedure STBIV;
begin
Hk[1]:=0.000000000000000000;
Hk[2]:=0.056262560526922147;
Hk[3]:=0.180240691736892365;
Hk[4]:=0.352624717113169337;
Hk[5]:=0.547153626330555383;
Hk[6]:=0.734210177215410532;
Hk[7]:=0.885320946839095768;
Hk[8]:=0.977520613561287501;
end;

Procedure BIV;
var Pi:longint;
begin
for i:=1 to Step do Thk[i]:=H*Hk[i];
for i:=1 to Step do
for Pi:=1 to Step do begin
if (Pi>1) and (i>Pi) then Ca[i,Pi]:=Ca[i-1,Pi-1]-Thk[i]*Ca[i-1,Pi];
if (Pi=1) and (i>1) then Ca[i,Pi]:=-Thk[i]*Ca[i-1,1];
if i=Pi then Ca[i,Pi]:=1; end;
end;


Procedure EVR;
begin
STBIV;
Np:=1;
repeat
for i:=1 to N do Y0[i]:=YY[i];
if h<0 then if t+h<TPR[NP] then h:=TPR[NP]-T;
if h>0 then if t+h>TPR[NP] then h:=TPR[NP]-T; t:=t+h; BIV;
DIFFUN1; for i:=1 to N do begin D02[i]:=Y0[i]+D[i]*h; D0[i]:=D[i]; end;
XYHK(Thk[2]); DIFFUN1; AMA(1); AVA(1);
XYHK(Thk[3]); DIFFUN1; AMA(2); AVA(2);
XYHK(Thk[4]); DIFFUN1; AMA(3); AVA(3);
XYHK(Thk[5]); DIFFUN1; AMA(4); AVA(4);
XYHK(Thk[6]); DIFFUN1; AMA(5); AVA(5);
XYHK(Thk[7]); DIFFUN1; AMA(6); AVA(6);
XYHK(Thk[8]); DIFFUN1; AMA(7); AVA(7);
XYHK(H);
{ for i:=1 to N do Epsm:=Epsm+(abs(Yy[i]-D02[i]));
if Epsm>Eps then begin t:=t-h; h:=0.2*h; pp:=pp-1; for i:=1 to N do Yy[i]:=Y0[i]; end;
if 0.5*Epsm<Eps then H:=2*h;}
if T=TPR[NP] then begin OUTPR; if NP<>NPR then NP:=NP+1; end;
PP:=PP+1; E_Kinet; U_poten; if round(pp/hi)=pp/hi then GRG1;{Графика}
until T=TPR[NPR];
end;