implementation
//PROGRAM D12R2
//Driver for routine RKDUMB
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;  NSTEP = 150;
var
  F:TextFile;  Y:matrx2;
  VSTART:array[0..NVAR] of real;  XX:array[0..200] of real;
  I,J:integer;   X1,X2:real;
begin
  SetLength(Y,11,201);
  X1:=1;
  VSTART[1]:= BESSJ0(X1);
  VSTART[2]:= BESSJ1(X1);
  VSTART[3]:= BESSJ(2, X1);
  VSTART[4]:= BESSJ(3, X1);
  X2:= 20;
  RKDUMB(VSTART, NVAR, X1, X2, NSTEP, XX, Y);
  //輸出計算結果到文件
  AssignFile(F, 'd:\delphi_shu\p12\d12r2.dat');
  Rewrite(F);
  Writeln(F);
  Writeln(F,'      X        Integrated      BESSJ3');
  For I:= 1 To Round(NSTEP / 10) do
  begin
    J:= 10 * I;
    Writeln(F,Format(s2,[XX[J]]),Format(s1,[Y[4, J]]),
                                    Format(s1,[BESSJ(3, XX[J])]));
  end;
  CloseFile(F);
  //屏幕顯示計算結果
  memo1.Lines.LoadFromFile('d:\delphi_shu\p12\d12r2.dat');
end;