implementation
//PROGRAM D12R5
//Driver for routine MMID
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 = '%4.2f';  NVAR = 4;  X1 = 1.0;  HTOT = 0.5;
var
  F:TextFile;   Y,YOUT,DYDX:array[0..NVAR] of real;
  I,II:integer;   XF,B1,B2,B3,B4:real;
begin
  Y[1]:= BESSJ0(X1);
  Y[2]:= BESSJ1(X1);
  Y[3]:= BESSJ(2, X1);
  Y[4]:= BESSJ(3, X1);
  DYDX[1]:= -Y[2];
  DYDX[2]:= Y[1] - Y[2];
  DYDX[3]:= Y[2] - 2*Y[3];
  DYDX[4]:= Y[3] - 3*Y[4];
  XF:= X1 + HTOT;
  B1:= BESSJ0(XF);
  B2:= BESSJ1(XF);
  B3:= BESSJ(2, XF);
  B4:= BESSJ(3, XF);
  //輸出計算結果到文件
  AssignFile(F, 'd:\delphi_shu\p12\d12r5.dat');
  Rewrite(F);
  Writeln(F);
  Writeln(F,'First four Bessl functions');
  For II:= 0 To (50-5) div 5 do
  begin
    I:= 5 + II*5;
    MMID(Y, DYDX, NVAR, X1, HTOT, I, YOUT);
    Writeln(F);
    Writeln(F,'X= ', Format(s2,[X1]),' to ',
       Format(s2,[X1 + HTOT]),' in ', FormatFloat('#0',I), ' steps');
    Writeln(F,'     Integration      BESSJ');
    Writeln(F, Format(s1,[YOUT[1]]), Format(s1,[B1]));
    Writeln(F, Format(s1,[YOUT[2]]), Format(s1,[B2]));
    Writeln(F, Format(s1,[YOUT[3]]), Format(s1,[B3]));
    Writeln(F, Format(s1,[YOUT[4]]), Format(s1,[B4]));
  end;
  CloseFile(F);
  //屏幕顯示計算結果
  memo1.Lines.LoadFromFile('d:\delphi_shu\p12\d12r5.dat');
end;