implementation
//PROGRAM D12R4
//Driver for routine ODEINT
uses
  unit2;
  {$R *.DFM}
procedure DERIVS(X:real; Y:array of real;var DYDX:array of real);
begin
    DYDX[1]:= -Y[2];
    DYDX[2]:= Y[1] - (1 / X) * Y[2];
    DYDX[3]:= Y[2] - (2 / X) * Y[3];
    DYDX[4]:= Y[3] - (3 / X) * Y[4];
end;
procedure TForm1.Button1Click(Sender: TObject);
const
  s1='%14.6f';  s2='%10.4f';NVAR = 4;
var
  F:TextFile;  YP:matrx2;
  VSTART:array[0..NVAR] of real;  XP:array[0..200] of real;
  I,NOK,NBAD:integer;   X1,X2,H1,EPS,HMIN:real;
begin
  SetLength(YP,11,201);
  X1:=1.0;
  X2:=10.0;
  VSTART[1]:= BESSJ0(X1);
  VSTART[2]:= BESSJ1(X1);
  VSTART[3]:= BESSJ(2, X1);
  VSTART[4]:= BESSJ(3, X1);
  EPS:=0.0001;
  H1:= 0.1;
  HMIN:= 0;
  KMAX:= 100;
  DXSAV:= (X2 - X1) / 20;
  ODEINT(VSTART, NVAR, X1, X2, EPS, H1, HMIN, NOK, NBAD, XP, YP);
  //輸出計算結(jié)果到文件
  AssignFile(F, 'd:\delphi_shu\p12\d12r4.dat');
  Rewrite(F);
  Writeln(F);
  Writeln(F,'Successful step:             ',FormatFloat('##0',NOK));
  Writeln(F,'Bad step:                    ',FormatFloat('##0',NBAD));
  Writeln(F,'Stored intermediate values:  ',FormatFloat('##0',KOUNT));
  Writeln(F);
  Writeln(F,'      X         Integral      BESSJ(3,X)');
  For I:= 1 To KOUNT do
        Writeln(F,Format(s2,[XP[I]]),Format(s1,[YP[4, I]]),
                                 Format(s1,[BESSJ(3, XP[I])]));
  CloseFile(F);
  //屏幕顯示計算結(jié)果
  memo1.Lines.LoadFromFile('d:\delphi_shu\p12\d12r4.dat');
end;