implementation
//PROGRAM D12R1
//Driver for routine RK4
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='%12.6f';  N = 4;
var
  F:TextFile;
  Y,DYDX,YOUT:array[0..N] of real;  I,J:integer;   X,H:real;
begin
  X:=1;
  Y[1]:=BESSJ0(X);
  Y[2]:=BESSJ1(X);
  Y[3]:=BESSJ(2, X);
  Y[4]:=BESSJ(3, X);
  DYDX[1]:=-Y[2];
  DYDX[2]:=Y[1] - Y[2];
  DYDX[3]:=Y[2] - 2* Y[3];
  DYDX[4]:=Y[3] - 3* Y[4];
  //輸出計算結果到文件
  AssignFile(F, 'd:\delphi_shu\p12\d12r1.dat');
  Rewrite(F);
  Writeln(F);
  Writeln(F,'Bessel Function:     J0          J1          J2         J3');
  For I:=1 To 5 do
  begin
    H:=0.2 * I;
    RK4(Y, DYDX, N, X, H, YOUT);
    Writeln(F);
    Writeln(F,'For a step size of: ', Format('%4.2f',[H]));
    Writeln(F,'      RK4:    ',Format(s1,[YOUT[1]]),Format(s1,[YOUT[2]]),
               Format(s1,[YOUT[3]]),Format(s1,[YOUT[4]]));
    Writeln(F,'      Actual: ',Format(s1,[BESSJ0(X + H)]),
               Format(s1,[BESSJ1(X + H)]),Format(s1,[BESSJ(2, X + H)]),
               Format(s1,[BESSJ(3, X + H)]));
  end;
  CloseFile(F);
  //屏幕顯示計算結果
  memo1.Lines.LoadFromFile('d:\delphi_shu\p12\d12r1.dat');
end;