program IntervalMultistepMethods;
{$APPTYPE CONSOLE}
uses System.SysUtils, Winapi.Windows, IntervalArithmetic32and64;

type interval_function = function (const T, Y : interval) : interval;

procedure AB1 (const T0, Y0, Fdelta : interval;
               const h              : Extended;
               const FTY, PSI       : interval_function;
               out T1, Y1           : interval);
// explicit interval one-step method of Adams-Bashforth type
var i1, i2, ih, hh, T, Y : interval;
    s                    : string;
begin
  Str (h:26, s);
  ih.a:=left_read(s);
  ih.b:=right_read(s);
  hh.a:=0;
  hh.b:=ih.b;
  T1:=T0+ih;
  T:=T0+hh;
  Y:=Y0+hh*Fdelta;
  i1:=PSI(T,Y);
  i1:=ih*ih*i1;
  i2:=FTY(T0,Y0);
  i2:=ih*i2;
  i2:=i2+i1/2;
  Y1:=Y0+i2
end {AB1};

procedure N1 (const T0, Y0, Y1, Fdelta : interval;
              const h                  : Extended;
              const FTY, PSI           : interval_function;
              out T1, T2, Y2           : interval);
// explicit interval two-step method of Nystrom type
var i1, i2, ih, hh, T, Y : interval;
    s                    : string;
begin
  Str (h:26, s);
  ih.a:=left_read(s);
  ih.b:=right_read(s);
  hh.a:=0;
  hh.b:=ih.b;
  T1:=T0+ih;
  T2:=T1+ih;
  T:=T1+hh;
  Y:=Y1+hh*Fdelta;
  i1:=PSI(T,Y);
  i1:=ih*ih*(i1-i1);
  i2:=FTY(T1,Y1);
  i2:=2*ih*i2;
  i2:=i2+i1/2;
  Y2:=Y0+i2
end {N1};

procedure AB2 (const T0, Y0, Y1, Fdelta : interval;
               const h                  : Extended;
               const FTY, PSI           : interval_function;
               out T1, T2, Y2           : interval);
// explicit interval two-step method of Adams-Bashforth type
var i1, i2, i3, ih, hh, T, Y : interval;
    s                        : string;
begin
  Str (h:26, s);
  ih.a:=left_read(s);
  ih.b:=right_read(s);
  Str (-h:26, s);
  hh.a:=left_read(s);
  hh.b:=ih.b;
  T1:=T0+ih;
  T2:=T1+ih;
  T:=T1+hh;
  Y:=Y1+hh*Fdelta;
  i1:=PSI(T,Y);
  i2:=ih*ih*ih;
  i2:=5*i2*i1;
  i1:=i2/12;
  i2:=FTY(T1,Y1);
  i3:=FTY(T0,Y0);
  i2:=3*i2-i3;
  i2:=ih*i2/2;
  Y2:=Y1+(i2+i1)
end {AB2};

procedure N2 (const T0, Y0, Y1, Fdelta : interval;
              const h                  : Extended;
              const FTY, PSI           : interval_function;
              out T1, T2, Y2           : interval);
//  explicit interval two-step method of Nystrom type
//  (in conventional case the same as for k=1)
var i1, i2, ih, hh, T, Y : interval;
    s                    : string;
begin
  Str (h:26, s);
  ih.a:=left_read(s);
  ih.b:=right_read(s);
  Str (-h:26, s);
  hh.a:=left_read(s);
  hh.b:=ih.b;
  T1:=T0+ih;
  T2:=T1+ih;
  T:=T1+hh;
  Y:=Y1+hh*Fdelta;
  i1:=PSI(T,Y);
  i1:=5*i1-i1;
  i2:=ih*ih*ih;
  i2:=i2*i1/12;
  i1:=FTY(T1,Y1);
  i1:=2*ih*i1;
  Y2:=Y0+(i1+i2)
end {N2};

procedure AB3 (const T0, Y0, Y1, Y2, Fdelta : interval;
               const h                      : Extended;
               const FTY, PSI               : interval_function;
               out T1, T2, T3, Y3           : interval);
// explicit interval three-step method of Adams_Bashforth type
var i1, i2, i3, i4, ih, hh, T, Y : interval;
    s                            : string;
begin
  Str (h:26, s);
  ih.a:=left_read(s);
  ih.b:=right_read(s);
  Str (-2*h:26, s);
  hh.a:=left_read(s);
  hh.b:=ih.b;
  T1:=T0+ih;
  T2:=T1+ih;
  T3:=T2+ih;
  T:=T2+hh;
  Y:=Y2+hh*Fdelta;
  i1:=PSI(T,Y);
  i2:=ih*ih;
  i2:=i2*i2;
  i2:=3*i2*i1;
  i1:=i2/8;
  i2:=FTY(T2,Y2);
  i3:=FTY(T1,Y1);
  i4:=FTY(T0,Y0);
  i2:=23*i2-16*i3;
  i2:=i2+5*i4;
  i2:=ih*i2/12;
  Y3:=Y2+(i2+i1)
end {AB3};

procedure N3 (const T0, Y0, Y1, Y2, Fdelta : interval;
              const h                      : Extended;
              const FTY, PSI               : interval_function;
              out T1, T2, T3, Y3           : interval);
// explicit interval three-step method of Nystrom type
var i1, i2, i3, i4, ih, hh, T, Y : interval;
    s                            : string;
begin
  Str (h:26, s);
  ih.a:=left_read(s);
  ih.b:=right_read(s);
  Str (-2*h:26, s);
  hh.a:=left_read(s);
  hh.b:=ih.b;
  T1:=T0+ih;
  T2:=T1+ih;
  T3:=T2+ih;
  T:=T2+hh;
  Y:=Y2+hh*Fdelta;
  i1:=PSI(T,Y);
  i1:=9*i1-i1;
  i2:=ih*ih;
  i2:=i2*i2;
  i1:=i2*i1/24;
  i2:=FTY(T2,Y2);
  i3:=FTY(T1,Y1);
  i4:=FTY(T0,Y0);
  i2:=7*i2-2*i3;
  i2:=i2+i4;
  i2:=ih*i2/3;
  Y3:=Y1+(i2+i1)
end {N3};

procedure AB4 (const T0, Y0, Y1, Y2, Y3, Fdelta : interval;
               const h                          : Extended;
               const FTY, PSI                   : interval_function;
               out T1, T2, T3, T4, Y4           : interval);
// explicit interval four-step method of Adams-Bashforth type
var i1, i2, i3, i4, i5, ih, hh, T, Y : interval;
    s                                : string;
begin
  Str (h:26, s);
  ih.a:=left_read(s);
  ih.b:=right_read(s);
  Str (-3*h:26, s);
  hh.a:=left_read(s);
  hh.b:=ih.b;
  T1:=T0+ih;
  T2:=T1+ih;
  T3:=T2+ih;
  T4:=T3+ih;
  T:=T3+hh;
  Y:=Y3+hh*Fdelta;
  i1:=PSI(T,Y);
  i2:=ih*ih;
  i2:=i2*i2;
  i2:=i2*ih;
  i2:=271*i2*i1;
  i1:=i2/720;
  i2:=FTY(T3,Y3);
  i3:=FTY(T2,Y2);
  i4:=FTY(T1,Y1);
  i5:=FTY(T0,Y0);
  i2:=55*i2-59*i3;
  i2:=i2+37*i4;
  i2:=i2-9*i5;
  i2:=ih*i2/24;
  Y4:=Y3+(i2+i1)
end {AB4};

procedure N4 (const T0, Y0, Y1, Y2, Y3, Fdelta : interval;
              const h                          : Extended;
              const FTY, PSI                   : interval_function;
              out T1, T2, T3, T4, Y4           : interval);
// explicit interval four-step method of Nystrom type
var i1, i2, i3, i4, i5, ih, hh, T, Y : interval;
    s                                : string;
begin
  Str (h:26, s);
  ih.a:=left_read(s);
  ih.b:=right_read(s);
  Str (-3*h:26, s);
  hh.a:=left_read(s);
  hh.b:=ih.b;
  T1:=T0+ih;
  T2:=T1+ih;
  T3:=T2+ih;
  T4:=T3+ih;
  T:=T3+hh;
  Y:=Y3+hh*Fdelta;
  i1:=PSI(T,Y);
  i1:=251*i1-19*i1;
  i2:=ih*ih;
  i2:=i2*i2;
  i2:=i2*ih*i1;
  i1:=i2/720;
  i2:=FTY(T3,Y3);
  i3:=FTY(T2,Y2);
  i4:=FTY(T1,Y1);
  i5:=FTY(T0,Y0);
  i2:=8*i2-5*i3;
  i2:=i2+4*i4;
  i2:=i2-i5;
  i2:=ih*i2/3;
  Y4:=Y2+(i2+i1)
end {N4};

procedure AM1 (const T0, Y0, Fdelta : interval;
               const h              : Extended;
               const FTY, PSI       : interval_function;
               out find             : Boolean;
               out T1, Y1           : interval;
               out l                : Integer);
// implicit interval one-step method of Adams-Moulton type
var i1, i2, i3, ih, hh, T, Y, Y2 : interval;
    s                            : string;
begin
  Str (h:26, s);
  ih.a:=left_read(s);
  ih.b:=right_read(s);
  Str (-h:26, s);
  hh.a:=left_read(s);
  hh.b:=0;
  T1:=T0+ih;
  Y1:=Y0;
  T:=T1+hh;
//  iterations
  l:=0;  // number of iterations
  repeat
    l:=l+1;
    find:=True;
    Y:=Y1+hh*Fdelta;
    i1:=PSI(T,Y);
    i2:=ih*ih*ih;
    i1:=i2*i1/12;
    i2:=FTY(T1,Y1);
    i3:=FTY(T0,Y0);
    i2:=i2+i3;
    i2:=ih*i2/2;
    Y2:=Y0+(i2-i1);
    if (Y1.a=0) or (Y1.b=0)
      then if (Abs(Y2.a-Y1.a)>=1e-18) or (Abs(Y2.b-Y1.b)>=1e-18)
             then find:=False  // 1e-18 - accuracy
             else
      else if (Abs((Y2.a-Y1.a)/Y1.a)>=1e-18) or (Abs((Y2.b-Y1.b)/Y1.b)>=1e-18)
             then find:=False;
    Y1:=Y2;
  until (l=20) or find  // 20 - max number of iterations
end {AM1};

procedure MS1 (const T0, Y0, Y1, Fdelta : interval;
               const h                  : Extended;
               const FTY, PSI           : interval_function;
               out find                 : Boolean;
               out T1, T2, Y2           : interval;
               out l                    : Integer);
// implicit interval two-step method of Milne-Simpson type (k = 1)
var i1, i2, ih, hh, T, Y, Y3 : interval;
    s                        : string;
begin
  Str (h:26, s);
  ih.a:=left_read(s);
  ih.b:=right_read(s);
  Str (-h:26, s);
  hh.a:=left_read(s);
  hh.b:=0;
  T1:=T0+ih;
  T2:=T1+ih;
  Y2:=Y1;
  T:=T2+hh;
//  iterations
  l:=0;  // number of iterations
  repeat
    l:=l+1;
    find:=True;
    Y:=Y2+hh*Fdelta;
    i1:=PSI(T,Y);
    i1:=5*i1-i1;
    i2:=ih*ih*ih;
    i1:=i2*i1/12;
    i2:=FTY(T1,Y1);
    i2:=2*ih*i2;
    Y3:=Y0+(i2+i1);
    if (Y2.a=0) or (Y2.b=0)
      then if (Abs(Y3.a-Y2.a)>=1e-18) or (Abs(Y3.b-Y2.b)>=1e-18)
             then find:=False  // 1e-18 - accuracy
             else
      else if (Abs((Y3.a-Y2.a)/Y2.a)>=1e-18) or (Abs((Y3.b-Y2.b)/Y2.b)>=1e-18)
             then find:=False;
    Y2:=Y3;
  until (l=20) or find  // 20 - max number of iterations
end {MS1};

procedure AM2 (const T0, Y0, Y1, Fdelta : interval;
               const h                  : Extended;
               const FTY, PSI           : interval_function;
               out find                 : Boolean;
               out T1, T2, Y2           : interval;
               out l                    : Integer);
// implicit interval two-step method of Adams-Moulton type
var i1, i2, i3, i4, ih, hh, T, Y, Y3 : interval;
    s                                : string;
begin
  Str (h:26, s);
  ih.a:=left_read(s);
  ih.b:=right_read(s);
  Str (-2*h:26, s);
  hh.a:=left_read(s);
  hh.b:=0;
  T1:=T0+ih;
  T2:=T1+ih;
  Y2:=Y1;
  T:=T2+hh;
// iterations
  l:=0;  // number of iterations
  repeat
    l:=l+1;
    find:=True;
    Y:=Y2+hh*Fdelta;
    i1:=PSI(T,Y);
    i2:=ih*ih;
    i2:=i2*i2;
    i1:=i2*i1/24;
    i2:=FTY(T2,Y2);
    i3:=FTY(T1,Y1);
    i4:=FTY(T0,Y0);
    i2:=5*i2+8*i3;
    i2:=i2-i4;
    i2:=ih*i2/12;
    Y3:=Y1+(i2-i1);
    if (Y2.a=0) or (Y2.b=0)
      then if (Abs(Y3.a-Y2.a)>=1e-18) or (Abs(Y3.b-Y2.b)>=1e-18)
             then find:=False  // 1e-18 - accuracy
             else
      else if (Abs((Y3.a-Y2.a)/Y2.a)>=1e-18) or (Abs((Y3.b-Y2.b)/Y2.b)>=1e-18)
             then find:=False;
    Y2:=Y3;
  until (l=20) or find  // 20 - max number of iterations
end {AM2};

procedure MS2 (const T0, Y0, Y1, Fdelta : interval;
               const h                  : Extended;
               const FTY, PSI           : interval_function;
               out find                 : Boolean;
               out T1, T2, Y2           : interval;
               out l                    : Integer);
// implicit interval two-step method of Milne-Simpson type (k = 2)
var i1, i2, i3, ih, hh, T, Y, Y3 : interval;
    s                            : string;
begin
  Str (h:26, s);
  ih.a:=left_read(s);
  ih.b:=right_read(s);
  Str (-2*h:26, s);
  hh.a:=left_read(s);
  hh.b:=0;
  T1:=T0+ih;
  T2:=T1+ih;
  Y2:=Y1;
  T:=T2+hh;
// iterations
  l:=0;  // number of iterations
  repeat
    l:=l+1;
    find:=True;
    Y:=Y2+hh*Fdelta;
    i1:=PSI(T,Y);
    i1:=i1-i1;
    i2:=ih*ih;
    i2:=i2*i2;
    i1:=i2*i1/24;
    i2:=FTY(T2,Y2);
    i3:=FTY(T1,Y1);
    i2:=i2+4*i3;
    i3:=FTY(T0,Y0);
    i2:=i2+i3;
    i2:=ih*i2/3;
    Y3:=Y0+(i2+i1);
    if (Y2.a=0) or (Y2.b=0)
      then if (Abs(Y3.a-Y2.a)>=1e-18) or (Abs(Y3.b-Y2.b)>=1e-18)
             then find:=False  // 1e-18 - accuracy
             else
      else if (Abs((Y3.a-Y2.a)/Y2.a)>=1e-18) or (Abs((Y3.b-Y2.b)/Y2.b)>=1e-18)
             then find:=False;
    Y2:=Y3;
  until (l=20) or find  // 20 - max number of iterations
end {MS2};

procedure AM3 (const T0, Y0, Y1, Y2, Fdelta : interval;
               const h                      : Extended;
               const FTY, PSI               : interval_function;
               out find                     : Boolean;
               out T1, T2, T3, Y3           : interval;
               out l                        : Integer);
// implicit interval three-step method of Adams-Moulton type
var i1, i2, i3, i4, i5, ih, hh, T, Y, Y4 : interval;
    s                                    : string;
begin
  Str (h:26, s);
  ih.a:=left_read(s);
  ih.b:=right_read(s);
  Str (-3*h:26, s);
  hh.a:=left_read(s);
  hh.b:=0;
  T1:=T0+ih;
  T2:=T1+ih;
  T3:=T2+ih;
  Y3:=Y2;
  T:=T3+hh;
// iterations
  l:=0;  // number of iterations
  repeat
    l:=l+1;
    find:=True;
    Y:=Y3+hh*Fdelta;
    i1:=PSI(T,Y);
    i3:=ih*ih;
    i2:=i3*ih;
    i2:=i2*i3;
    i2:=19*i2*i1;
    i1:=i2/720;
    i2:=FTY(T3,Y3);
    i3:=FTY(T2,Y2);
    i4:=FTY(T1,Y1);
    i5:=FTY(T0,Y0);
    i2:=9*i2+19*i3;
    i2:=i2-5*i4;
    i2:=i2+i5;
    i2:=ih*i2/24;
    Y4:=Y2+(i2-i1);
    if (Y3.a=0) or (Y3.b=0)
      then if (Abs(Y4.a-Y3.a)>=1e-18) or (Abs(Y4.b-Y3.b)>=1e-18)
             then find:=False  // 1e-18 - accuracy
             else
      else if (Abs((Y4.a-Y3.a)/Y3.a)>=1e-18) or (Abs((Y4.b-Y3.b)/Y3.b)>=1e-18)
             then find:=False;
    Y3:=Y4;
  until (l=20) or find  // 20 - max number of iterations
end {AM3};

procedure MS3 (const T0, Y0, Y1, Fdelta : interval;
               const h                  : Extended;
               const FTY, PSI           : interval_function;
               out find                 : Boolean;
               out T1, T2, Y2           : interval;
               out l                    : Integer);
// implicit interval two-step method of Milne-Simpson type (k = 3)
// (in conventional case the same as for k=2)
var i1, i2, i3, ih, hh, T, Y, Y3 : interval;
    s                            : string;
begin
  Str (h:26, s);
  ih.a:=left_read(s);
  ih.b:=right_read(s);
  Str (-3*h:26, s);
  hh.a:=left_read(s);
  hh.b:=0;
  T1:=T0+ih;
  T2:=T1+ih;
  Y2:=Y1;
  T:=T2+hh;
// iterations
  l:=0;  // number of iterations
  repeat
    l:=l+1;
    find:=True;
    Y:=Y2+hh*Fdelta;
    i1:=PSI(T,Y);
    i1:=11*i1-19*i1;
    i2:=ih*ih;
    i2:=i2*i2;
    i2:=i2*ih;
    i1:=i2*i1/720;
    i2:=FTY(T2,Y2);
    i3:=FTY(T1,Y1);
    i2:=i2+4*i3;
    i3:=FTY(T0,Y0);
    i2:=i2+i3;
    i2:=ih*i2/3;
    Y3:=Y0+(i2+i1);
    if (Y2.a=0) or (Y2.b=0)
      then if (Abs(Y3.a-Y2.a)>=1e-18) or (Abs(Y3.b-Y2.b)>=1e-18)
             then find:=False  // 1e-18 - accuracy
             else
      else if (Abs((Y3.a-Y2.a)/Y2.a)>=1e-18) or (Abs((Y3.b-Y2.b)/Y2.b)>=1e-18)
             then find:=False;
    Y2:=Y3;
  until (l=20) or find  // 20 - max number of iterations
end {MS3};

procedure ABM1 (const T0, Y0, Fdelta  : interval;
                const h               : Extended;
                const FTY, PSI1, PSI2 : interval_function;
                out find              : Boolean;
                out T1, Y1            : interval;
                out l                 : Integer);
// implicit interval one-step predictor-corrector method of
// Adams-Bashforth-Moulton type
var i1, i2, i3, ih, hh, hh1, T, Y, Y2 : interval;
    s                                 : string;
begin
  Str (h:26, s);
  ih.a:=left_read(s);
  ih.b:=right_read(s);
  hh.a:=0;
  hh.b:=ih.b;
  Str (-h:26, s);
  hh1.a:=left_read(s);
  hh1.b:=0;
  T1:=T0+ih;
  T:=T0+hh;
  Y:=Y0+hh*Fdelta;
  i1:=PSI1(T,Y);
  i1:=ih*ih*i1;
  i2:=FTY(T0,Y0);
  i2:=ih*i2;
  i2:=i2+i1/2;
  Y1:=Y0+i2;
  T:=T1+hh1;
// iterations
  l:=0;  // number of iterations
  repeat
    l:=l+1;
    find:=True;
    Y:=Y1+hh1*Fdelta;
    i1:=PSI2(T,Y);
    i2:=ih*ih*ih;
    i1:=i2*i1/12;
    i2:=FTY(T1,Y1);
    i3:=FTY(T0,Y0);
    i2:=i2+i3;
    i2:=ih*i2/2;
    Y2:=Y0+(i2-i1);
    if (Y1.a=0) or (Y1.b=0)
      then if (Abs(Y2.a-Y1.a)>=1e-18) or (Abs(Y2.b-Y1.b)>=1e-18)
             then find:=False  // 1e-18 - accuracy
             else
      else if (Abs((Y2.a-Y1.a)/Y1.a)>=1e-18) or (Abs((Y2.b-Y1.b)/Y1.b)>=1e-18)
             then find:=False;
    Y1:=Y2;
  until (l=20) or find  // 20 - max number of iterations
end {ABM1};

procedure NMS1 (const T0, Y0, Y1, Fdelta  : interval;
                const h                   : Extended;
                const FTY, PSI1, PSI2     : interval_function;
                out find                  : Boolean;
                out T1, T2, Y2            : interval;
                out l                     : Integer);
// implicit interval two-step predictor-corrector method of
// Nystrom-Milne-Simpson type (k = 1)
var i1, i2, ih, hh, hh1, T, Y, Y3 : interval;
    s                             : string;
begin
  Str (h:26, s);
  ih.a:=left_read(s);
  ih.b:=right_read(s);
  hh.a:=0;
  hh.b:=ih.b;
  Str (-h:26, s);
  hh1.a:=left_read(s);
  hh1.b:=0;
  T1:=T0+ih;
  T2:=T1+ih;
  T:=T1+hh;
  Y:=Y1+hh*Fdelta;
  i1:=PSI1(T,Y);
  i1:=ih*ih*(i1-i1);
  i2:=FTY(T1,Y1);
  i2:=2*ih*i2;
  i2:=i2+i1/2;
  Y2:=Y0+i2;
  T:=T2+hh1;
// iterations
  l:=0;  // number of iterations
  repeat
    l:=l+1;
    find:=True;
    Y:=Y2+hh1*Fdelta;
    i1:=PSI2(T,Y);
    i1:=5*i1-i1;
    i2:=ih*ih*ih;
    i1:=i2*i1/12;
    i2:=FTY(T1,Y1);
    i2:=2*ih*i2;
    Y3:=Y0+(i2+i1);
    if (Y2.a=0) or (Y2.b=0)
      then if (Abs(Y3.a-Y2.a)>=1e-18) or (Abs(Y3.b-Y2.b)>=1e-18)
             then find:=False  // 1e-18 - accuracy
             else
      else if (Abs((Y3.a-Y2.a)/Y2.a)>=1e-18) or (Abs((Y3.b-Y2.b)/Y2.b)>=1e-18)
             then find:=False;
    Y2:=Y3
  until (l=20) or find  // 20 - max number of iterations
end {NM1};

procedure ABM2 (const T0, Y0, Y1, Fdelta  : interval;
                const h                   : Extended;
                const FTY, PSI1, PSI2     : interval_function;
                out find                  : Boolean;
                out T1, T2, Y2            : interval;
                out l                     : Integer);
// implicit interval two-step predictor-corrector method of
// Adams-Bashforth-Moulton type
var i1, i2, i3, i4, ih, hh, hh1, T, Y, Y3 : interval;
    s                                     : string;
begin
  Str (h:26, s);
  ih.a:=left_read(s);
  ih.b:=right_read(s);
  Str (-h:26, s);
  hh.a:=left_read(s);
  hh.b:=ih.b;
  Str (-2*h:26, s);
  hh1.a:=left_read(s);
  hh1.b:=0;
  T1:=T0+ih;
  T2:=T1+ih;
  T:=T1+hh;
  Y:=Y1+hh*Fdelta;
  i1:=PSI1(T,Y);
  i2:=ih*ih*ih;
  i2:=5*i2*i1;
  i1:=i2/12;
  i2:=FTY(T1,Y1);
  i3:=FTY(T0,Y0);
  i2:=3*i2-i3;
  i2:=ih*i2/2;
  Y2:=Y1+(i2+i1);
  T:=T2+hh1;
// iterations
  l:=0;  // number of iterations
  repeat
    l:=l+1;
    find:=True;
    Y:=Y2+hh1*Fdelta;
    i1:=PSI2(T,Y);
    i2:=ih*ih;
    i2:=i2*i2;
    i1:=i2*i1/24;
    i2:=FTY(T2,Y2);
    i3:=FTY(T1,Y1);
    i4:=FTY(T0,Y0);
    i2:=5*i2+8*i3;
    i2:=i2-i4;
    i2:=ih*i2/12;
    Y3:=Y1+(i2-i1);
    if (Y2.a=0) or (Y2.b=0)
      then if (Abs(Y3.a-Y2.a)>=1e-18) or (Abs(Y3.b-Y2.b)>=1e-18)
             then find:=False  // 1e-18 - accuracy
             else
      else if (Abs((Y3.a-Y2.a)/Y2.a)>=1e-18) or (Abs((Y3.b-Y2.b)/Y2.b)>=1e-18)
             then find:=False;
    Y2:=Y3;
  until (l=20) or find  // 20 - max number of iterations
end {ABM2};

procedure NMS2 (const T0, Y0, Y1, Fdelta  : interval;
                const h                   : Extended;
                const FTY, PSI1, PSI2     : interval_function;
                out find                  : Boolean;
                out T1, T2, Y2            : interval;
                out l                     : Integer);
// implicit interval two-step predictor-corrector method of
// Nystrom-Milne-Simpson type (k = 2)
var i1, i2, i3, ih, hh, hh1, T, Y, Y3 : interval;
    s                                 : string;
begin
  Str (h:26, s);
  ih.a:=left_read(s);
  ih.b:=right_read(s);
  Str (-h:26, s);
  hh.a:=left_read(s);
  hh.b:=ih.b;
  Str (-2*h:26, s);
  hh1.a:=left_read(s);
  hh1.b:=0;
  T1:=T0+ih;
  T2:=T1+ih;
  T:=T1+hh;
  Y:=Y1+hh*Fdelta;
  i1:=PSI1(T,Y);
  i1:=5*i1-i1;
  i2:=ih*ih*ih;
  i2:=i2*i1/12;
  i1:=FTY(T1,Y1);
  i1:=2*ih*i1;
  Y2:=Y0+(i1+i2);
  T:=T2+hh1;
// iterations
  l:=0;  // number of iterations
  repeat
    l:=l+1;
    find:=True;
    Y:=Y2+hh1*Fdelta;
    i1:=PSI2(T,Y);
    i1:=i1-i1;
    i2:=ih*ih;
    i2:=i2*i2;
    i1:=i2*i1/24;
    i2:=FTY(T2,Y2);
    i3:=FTY(T1,Y1);
    i2:=i2+4*i3;
    i3:=FTY(T0,Y0);
    i2:=i2+i3;
    i2:=ih*i2/3;
    Y3:=Y0+(i2+i1);
    if (Y2.a=0) or (Y2.b=0)
      then if (Abs(Y3.a-Y2.a)>=1e-18) or (Abs(Y3.b-Y2.b)>=1e-18)
             then find:=False  // 1e-18 - accuracy
             else
      else if (Abs((Y3.a-Y2.a)/Y2.a)>=1e-18) or (Abs((Y3.b-Y2.b)/Y2.b)>=1e-18)
             then find:=False;
    Y2:=Y3
  until (l=20) or find  // 20 - max number of iterations
end {NMS2};

procedure ABM3 (const T0, Y0, Y1, Y2, Fdelta  : interval;
                const h                       : Extended;
                const FTY, PSI1, PSI2         : interval_function;
                out find                      : Boolean;
                out T1, T2, T3, Y3            : interval;
                out l                         : Integer);
// implicit interval three-step predictor-corrector method of
// Adams-Bashforth-Moulton type
var i1, i2, i3, i4, i5, ih, hh, hh1, T, Y, Y4 : interval;
    s                                         : string;
begin
  Str (h:26, s);
  ih.a:=left_read(s);
  ih.b:=right_read(s);
  Str (-2*h:26, s);
  hh.a:=left_read(s);
  hh.b:=ih.b;
  Str (-3*h:26, s);
  hh1.a:=left_read(s);
  hh1.b:=0;
  T1:=T0+ih;
  T2:=T1+ih;
  T3:=T2+ih;
  T:=T2+hh;
  Y:=Y2+hh*Fdelta;
  i1:=PSI1(T,Y);
  i2:=ih*ih;
  i2:=i2*i2;
  i2:=3*i2*i1;
  i1:=i2/8;
  i2:=FTY(T2,Y2);
  i3:=FTY(T1,Y1);
  i4:=FTY(T0,Y0);
  i2:=23*i2-16*i3;
  i2:=i2+5*i4;
  i2:=ih*i2/12;
  Y3:=Y2+(i2+i1);
  T:=T3+hh1;
// iterations
  l:=0;  // number of iterations
  repeat
    l:=l+1;
    find:=True;
    Y:=Y3+hh1*Fdelta;
    i1:=PSI2(T,Y);
    i3:=ih*ih;
    i2:=i3*ih;
    i2:=i2*i3;
    i2:=19*i2*i1;
    i1:=i2/720;
    i2:=FTY(T3,Y3);
    i3:=FTY(T2,Y2);
    i4:=FTY(T1,Y1);
    i5:=FTY(T0,Y0);
    i2:=9*i2+19*i3;
    i2:=i2-5*i4;
    i2:=i2+i5;
    i2:=ih*i2/24;
    Y4:=Y2+(i2-i1);
    if (Y3.a=0) or (Y3.b=0)
      then if (Abs(Y4.a-Y3.a)>=1e-18) or (Abs(Y4.b-Y3.b)>=1e-18)
             then find:=False  // 1e-18 - accuracy
             else
      else if (Abs((Y4.a-Y3.a)/Y3.a)>=1e-18) or (Abs((Y4.b-Y3.b)/Y3.b)>=1e-18)
             then find:=False;
    Y3:=Y4;
  until (l=20) or find  // 20 - max number of iterations
end;

procedure NMS3 (const T0, Y0, Y1, Y2, Fdelta  : interval;
                const h                       : Extended;
                const FTY, PSI1, PSI2         : interval_function;
                out find                      : Boolean;
                out T1, T2, T3, Y3            : interval;
                out l                         : Integer);
// implicit interval three-step predictor-corrector method of
// Nystrom-Milne-Simpson type
var i1, i2, i3, i4, ih, hh, hh1, T, Y, Y4 : interval;
    s                                     : string;
begin
  Str (h:26, s);
  ih.a:=left_read(s);
  ih.b:=right_read(s);
  Str (-2*h:26, s);
  hh.a:=left_read(s);
  hh.b:=ih.b;
  Str (-3*h:26, s);
  hh1.a:=left_read(s);
  hh1.b:=0;
  T1:=T0+ih;
  T2:=T1+ih;
  T3:=T2+ih;
  T:=T2+hh;
  Y:=Y2+hh*Fdelta;
  i1:=PSI1(T,Y);
  i1:=9*i1-i1;
  i2:=ih*ih;
  i2:=i2*i2;
  i1:=i2*i1/24;
  i2:=FTY(T2,Y2);
  i3:=FTY(T1,Y1);
  i4:=FTY(T0,Y0);
  i2:=7*i2-2*i3;
  i2:=i2+i4;
  i2:=ih*i2/3;
  Y3:=Y1+(i2+i1);
  T:=T3+hh1;
// iterations
  l:=0;  // number of iterations
  repeat
    l:=l+1;
    find:=True;
    Y:=Y3+hh1*Fdelta;
    i1:=PSI2(T,Y);
    i1:=11*i1-19*i1;
    i2:=ih*ih;
    i2:=i2*i2;
    i2:=i2*ih;
    i1:=i2*i1/720;
    i2:=FTY(T3,Y3);
    i3:=FTY(T2,Y2);
    i2:=i2+4*i3;
    i3:=FTY(T1,Y1);
    i2:=i2+i3;
    i2:=ih*i2/3;
    Y4:=Y1+(i2+i1);
    if (Y3.a=0) or (Y3.b=0)
      then if (Abs(Y4.a-Y3.a)>=1e-18) or (Abs(Y4.b-Y3.b)>=1e-18)
             then find:=False  // 1e-18 - accuracy
             else
      else if (Abs((Y4.a-Y3.a)/Y3.a)>=1e-18) or (Abs((Y4.b-Y3.b)/Y3.b)>=1e-18)
             then find:=False;
    Y3:=Y4;
  until (l=20) or find  // 20 - max number of iterations
end {NMS3};

var l, n, m, method, output_step      : Integer;
    amh, real_h, w                    : Extended;
    delta_t, delta_y, Fdelta, T0, Y0,
    T1, Y1, T2, Y2, T3, Y3, T4, Y4    : Interval;
    left, lib_name, right, strg       : string;
    s                                 : AnsiString;
    find                              : Boolean;
    z                                 : Char;
    plik                              : TextFile;
    ident                             : HMODULE;
    FTY, PSI1, PSI2                   : interval_function;

procedure Results (const T, Y : interval;
                   const kind : Char);
var hreal       : Extended;
    hrealstring : string;
begin
  hreal:=(T.a+T.b)/2;
  Str (hreal:3:1, hrealstring);
  iends_to_strings (T, left, right);
  if z='f'
    then Writeln (plik, 'T('+hrealstring+') = ['+left+','+right+']')
    else Writeln ('T('+hrealstring+') = ['+left+','+right+']');
  iends_to_strings (Y, left, right);
  if z='f'
    then Writeln (plik, 'Y('+hrealstring+') = ['+left+','+right+']')
    else Writeln ('Y('+hrealstring+') = ['+left+','+right+']');
  w:=int_width(Y);
  if z='f'
    then if kind='i'
           then Writeln (plik, 'width = ', w:11, ', number of iterations : ', l)
           else Writeln (plik, 'width = ', w:11)
    else if kind='i'
           then Writeln ('width = ', w:11, ', number of iterations : ', l)
           else Writeln ('width = ', w:11)
end {Results};

begin
  Writeln;
  Writeln ('Solving the initial value problem y'' = f(t,y(t)), y(0) = y0');
  Writeln;
  Writeln ('for delta_t = {t in R : 0 <= t <= a}');
  Writeln ('and delta_y = {y in R : b_left <= y <= b_right}');
  Writeln;
  Writeln ('by an interval multistep method.');
  Writeln;
  Writeln ('Choose an interval method (press Enter):');
  Readln;
  Writeln;
  Writeln (' 1 - explicit one-step first order of Adams-Bashforth type (k=1)');
  Writeln (' 2 - explicit two-step first order of Nystrom type (k=1)');
  Writeln (' 3 - explicit two-step second order of Adams-Bashforth type');
  Writeln (' 4 - explicit two-step second order of Nystrom type (k=2,');
  Writeln ('     in conventional case the same as for k=1)');
  Writeln (' 5 - explicit three-step third order of Adams-Bashforth type '
           +'(k=3)');
  Writeln (' 6 - explicit three-step third order of Nystrom type (k=3)');
  Writeln (' 7 - explicit four-step fourth order of Adams-Bashforth type '
           +'(k=4)');
  Writeln (' 8 - explicit four-step fourth order of Nystrom type (k=4)');
  Writeln (' 9 - implicit one-step second order of Adams-Moulton type (k=1)');
  Writeln ('10 - implicit two-step second order of Milne-Simpson type (k=1)');
  Writeln ('11 - implicit two-step third order of Adams-Moulton type (k=2)');
  Writeln ('12 - implicit two-step third order of Milne-Simpson type (k=2)');
  Writeln ('13 - implicit three-step fourth order of Adams-Moulton type (k=3)');
  Writeln ('14 - implicit two-step fourth order of Milne-Simpson type (k=3,');
  Writeln ('     in conventional case the same as for k=2)');
  Writeln ('15 - one-step predictor-corrector of Adams-Bashforth-Moulton type '
           +'(k=1)');
  Writeln ('16 - two-step predictor-corrector of Nystrom-Milne-Simpson type '
           +'(k=1)');
  Writeln ('17 - two-step predictor-corrector of Adams-Bashforth-Moulton type '
           +'(k=2)');
  Writeln ('18 - two-step predictor-corrector of Nystrom-Milne-Simpson type '
           +'(k=2)');
  Writeln ('19 - three-step predictor-corrector of Adams-Bashforth-Moulton '
           +'type (k=3)');
  Writeln ('20 - three-step predictor-corrector of Nystrom-Milne-Simpson type '
           +'(k=4)');
  Writeln;
  repeat
    Write ('Your choice (1 <= method <= 20) method = ');
    Readln (method)
  until (method>=1) and (method<=20);
  Writeln;
  Writeln ('You have choosen');
  Write ('the ');
  case method of
     1 : Write ('explicit one-step first order method of Adams-Bashforth');
     2 : Write ('explicit two-step first order method of Nystrom');
     3 : Write ('explicit two-step second order method of Adams-Bashforth');
     4 : Write ('explicit two-step second order method of Nystrom');
     5 : Write ('explicit three-step third order method of Adams-Bashforth');
     6 : Write ('explicit three-step third order method of Nystrom');
     7 : Write ('explicit four-step fourth order method of Adams-Bashforth');
     8 : Write ('explicit four-step fourth order method of Nystrom');
     9 : Write ('implicit one-step second order method of Adams-Moulton');
    10 : Write ('implicit two-step second order method of Milne-Simpson');
    11 : Write ('implicit two-step third order method of Adams-Moulton');
    12 : Write ('implicit two-step third order method of Milne-Simpson');
    13 : Write ('implicit three-step fourth order method of Adams-Moulton');
    14 : Write ('implicit two-step fourth order method of Milne-Simpson');
    15 : Write ('one-step predictor-corrector method of Adams-Bashforth-'
                +'Moulton');
    16 : Write ('two-step predictor-corrector method of Nystrom-Milne-Simpson');
    17 : Write ('two-step predictor-corrector method of Adams-Bashforth-'
                +'Moulton');
    18 : Write ('two-step predictor-corrector method of Nystrom-Milne-Simpson');
    19 : Write ('three-step predictor-corrector method of Adams-Bashforth-'
                +'Moulton');
    20 : Write ('three-step predictor-corrector method of Nystrom-Milne-'
                +'Simpson')
  end;
  Writeln (' type.');
  Writeln;
  Writeln ('This method requires a dynamic link library (DLL) with '
           +'definitions');
  Writeln ('of the following interval functions:');
  Writeln;
  Writeln ('F(T,Y)    - an interval extension of f(t,y(t)) with the header');
  Writeln ('            function FTY (const T, Y : interval) : interval;');
  Writeln;
  case method of
     1, 2  : begin
               Writeln ('PSI(T,Y)  - an interval extension of psi(t,y(t)) = '
                        +'f^(1)(t,y(t)) = y^(2)(t)');
               Writeln ('            with the header');
               Writeln ('            function Psi1Explicit (const T, Y : '
                        +'interval) : interval;')
             end;
     3, 4  : begin
               Writeln ('PSI(T,Y)  - an interval extension of psi(t,y(t)) = '
                        +'f^(2)(t,y(t)) = y^(3)(t)');
               Writeln ('            with the header');
               Writeln ('            function Psi2Explicit (const T, Y : '
                        +'interval) : interval;')
             end;
     5, 6  : begin
               Writeln ('PSI(T,Y)  - an interval extension of psi(t,y(t)) = '
                        +'f^(3)(t,y(t)) = y^(4)(t)');
               Writeln ('            with the header');
               Writeln ('            function Psi3Explicit (const T, Y : '
                        +'interval) : interval;')
             end;
     7, 8  : begin
               Writeln ('PSI(T,Y)  - an interval extension of psi(t,y(t)) = '
                        +'f^(4)(t,y(t)) = y^(5)(t)');
               Writeln ('            with the header');
               Writeln ('            function Psi4Explicit (const T, Y : '
                        +'interval) : interval;')
             end;
     9, 10 : begin
               Writeln ('PSI(T,Y)  - an interval extension of psi(t,y(t)) = '
                        +'f^(2)(t,y(t)) = y^(3)(t)');
               Writeln ('            with the header');
               Writeln ('            function Psi1Implicit (const T, Y : '
                        +'interval) : interval;')
             end;
    11, 12 : begin
               Writeln ('PSI(T,Y)  - an interval extension of psi(t,y(t)) = '
                        +'f^(3)(t,y(t)) = y^(4)(t)');
               Writeln ('            with the header');
               Writeln ('            function Psi2Implicit (const T, Y : '
                        +'interval) : interval;')
             end;
    13, 14 : begin
               Writeln ('PSI(T,Y)  - an interval extension of psi(t,y(t)) = '
                        +'f^(4)(t,y(t)) = y^(5)(t)');
               Writeln ('            with the header');
               Writeln ('            function Psi3Implicit (const T, Y : '
                        +'interval) : interval;')
             end;
    15, 16 : begin
               Writeln ('PSI1(T,Y) - an interval extension of psi(t,y(t)) = '
                        +'f^(1)(t,y(t)) = y^(2)(t)');
               Writeln ('            with the header');
               Writeln ('            function Psi1Explicit (const T, Y : '
                        +'interval) : interval;');
               Writeln;
               Writeln ('PSI2(T,Y) - an interval extension of psi(t,y(t)) = '
                        +'f^(2)(t,y(t)) = y^(3)(t)');
               Writeln ('            with the header');
               Writeln ('            function Psi1Implicit (const T, Y : '
                        +'interval) : interval;')
             end;
    17, 18 : begin
               Writeln ('PSI1(T,Y) - an interval extension of psi(t,y(t)) = '
                        +'f^(2)(t,y(t)) = y^(3)(t)');
               Writeln ('            with the header');
               Writeln ('            function Psi2Explicit (const T, Y : '
                        +'interval) : interval;');
               Writeln;
               Writeln ('PSI2(T,Y) - an interval extension of psi(t,y(t)) = '
                        +'f^(3)(t,y(t)) = y^(4)(t)');
               Writeln ('            with the header');
               Writeln ('            function Psi2Implicit (const T, Y : '
                        +'interval) : interval;')
             end;
    19, 20 : begin
               Writeln ('PSI1(T,Y) - an interval extension of psi(t,y(t)) = '
                        +'f^(3)(t,y(t)) = y^(4)(t)');
               Writeln ('            with the header');
               Writeln ('            function Psi3Explicit (const T, Y : '
                        +'interval) : interval;');
               Writeln;
               Writeln ('PSI2(T,Y) - an interval extension of psi(t,y(t)) = '
                        +'f^(4)(t,y(t)) = y^(5)(t)');
               Writeln ('            with the header');
               Writeln ('            function Psi3Implicit (const T, Y : '
                        +'interval) : interval;')
             end
  end;
  Writeln;
  repeat
    Write ('Do you have such a DLL (y/n)?: ');
    Readln (z)
  until (z='y') or (z='n');
  if z='n'
    then Exit;
  Writeln;
  repeat
    Writeln ('Enter the name of your DLL');
    Write ('(the extension .DLL will be added automatically): ');
    Readln (lib_name)
  until (lib_name<>'');
  lib_name:=lib_name+'.DLL';
  if not FileExists(lib_name)
    then begin
           Writeln;
           Writeln ('The file '+lib_name+' not found. The program will be '
                    +'closed.');
           Readln;
           Exit
         end;
  ident:=LoadLibrary(PWideChar(lib_name));
  if ident=0
    then begin
           Writeln;
           Writeln ('Your DLL has not been loaded. The program will be '
                    +'closed.');
           Readln;
           Exit
         end
    else begin
           Writeln;
           Writeln ('Your DLL has been loaded.')
         end;
  Writeln;
  find:=True;
  @FTY:=GetProcAddress(ident, 'FTY');
  if @FTY=nil
    then begin
           Writeln ('The function FTY has not been found in your DLL.');
           find:=False
         end
    else Writeln ('The function FTY has been found in your DLL.');
  if find
    then case method of
           1, 2  : begin
                     @PSI1:=GetProcAddress(ident,'Psi1Explicit');
                     if @PSI1=nil
                       then begin
                              Writeln ('The function Psi1Explicit has not been '
                                       +'found in your DLL.');
                              find:=False
                            end
                       else Writeln ('The function Psi1Explicit has been found '
                                     +'in your DLL.')
                   end;
           3, 4  : begin
                     @PSI1:=GetProcAddress(ident,'Psi2Explicit');
                     if @PSI1=nil
                       then begin
                              Writeln ('The function Psi2Explicit has not been '
                                       +'found in your DLL.');
                              find:=False
                            end
                       else Writeln ('The function Psi2Explicit has been found '
                                     +'in your DLL.')
                   end;
           5, 6  : begin
                     @PSI1:=GetProcAddress(ident,'Psi3Explicit');
                     if @PSI1=nil
                       then begin
                              Writeln ('The function Psi3Explicit has not been '
                                       +'found in your DLL.');
                              find:=False
                            end
                       else Writeln ('The function Psi3Explicit has been found '
                                     +'in your DLL.')
                   end;
           7, 8  : begin
                     @PSI1:=GetProcAddress(ident,'Psi4Explicit');
                     if @PSI1=nil
                       then begin
                              Writeln ('The function Psi4Explicit has not been '
                                       +'found in your DLL.');
                              find:=False
                            end
                       else Writeln ('The function Psi4Explicit has been found '
                                     +'in your DLL.')
                   end;
           9, 10 : begin
                     @PSI2:=GetProcAddress(ident,'Psi1Implicit');
                     if @PSI1=nil
                       then begin
                              Writeln ('The function Psi1Implicit has not been '
                                       +'found in your DLL.');
                              find:=False
                            end
                       else Writeln ('The function Psi1Implicit has been found '
                                     +'in your DLL.')
                   end;
          11, 12 : begin
                     @PSI2:=GetProcAddress(ident,'Psi2Implicit');
                     if @PSI1=nil
                       then begin
                              Writeln ('The function Psi2Implicit has not been '
                                       +'found in your DLL.');
                              find:=False
                            end
                       else Writeln ('The function Psi2Implicit has been found '
                                     +'in your DLL.')
                   end;
          13, 14 : begin
                     @PSI2:=GetProcAddress(ident,'Psi3Implicit');
                     if @PSI1=nil
                       then begin
                              Writeln ('The function Psi3Implicit has not been '
                                       +'found in your DLL.');
                              find:=False
                            end
                       else Writeln ('The function Psi3Implicit has been found '
                                     +'in your DLL.')
                   end;
          15, 16 : begin
                     @PSI1:=GetProcAddress(ident,'Psi1Explicit');
                     if @PSI1=nil
                       then begin
                              Writeln ('The function Psi1Explicit has not been '
                                       +'found in your DLL.');
                              find:=False
                            end
                       else begin
                              @PSI2:=GetProcAddress(ident,'Psi1Implicit');
                              if @PSI2=nil
                                then begin
                                       Writeln ('The function Psi1Implicit has '
                                                +'not been found in your DLL.');
                                       find:=False
                                     end
                                else Writeln ('The functions Psi1Explicit and '
                                              +'Psi1Implicit have been found '
                                              +'in your DLL.')
                            end
                   end;
          17, 18 : begin
                     @PSI1:=GetProcAddress(ident,'Psi2Explicit');
                     if @PSI1=nil
                       then begin
                              Writeln ('The function Psi2Explicit has not been '
                                       +'found in your DLL.');
                              find:=False
                            end
                       else begin
                              @PSI2:=GetProcAddress(ident,'Psi2Implicit');
                              if @PSI2=nil
                                then begin
                                       Writeln ('The function Psi2Implicit has '
                                                +'not been found in your DLL.');
                                       find:=False
                                     end
                                else Writeln ('The functions Psi2Explicit and '
                                              +'Psi2Implicit have been found '
                                              +'in your DLL.')
                            end
                   end;
          19, 20 : begin
                     @PSI1:=GetProcAddress(ident,'Psi3Explicit');
                     if @PSI1=nil
                       then begin
                              Writeln ('The function Psi3Explicit has not been '
                                       +'found in your DLL.');
                              find:=False
                            end
                       else begin
                              @PSI2:=GetProcAddress(ident,'Psi3Implicit');
                              if @PSI2=nil
                                then begin
                                       Writeln ('The function Psi3Implicit has '
                                                +'not been found in your DLL.');
                                       find:=False
                                     end
                                else Writeln ('The functions Psi3Explicit and '
                                              +'Psi3Implicit have been found '
                                              +'in your DLL.')
                            end
                   end
         end;
  if not find
    then begin
           Writeln ('The program will be closed.');
           FreeLibrary (ident);
           Readln;
           Exit
         end;
  Writeln;
  Writeln ('Enter data for your problem:');
  Writeln;
  Write ('delta_t = {t in R: 0 <= t <= a}                   a = ');
  Readln (w);
  amh:=w;
  Str (w:26, s);
  delta_t.a:=0;
  delta_t.b:=left_read(s);
  Write ('delta_y = (y in R: b_left <= y <= b_right)   b_left = ');
  Readln (w);
  Str (w:26, s);
  delta_y.a:=left_read(s);
  Write ('                                            b_right = ');
  Readln (w);
  Str (w:26, s);
  delta_y.b:=right_read(s);
  Fdelta:=FTY(delta_t,delta_y);
  Write ('step size                                         h = ');
  Readln (real_h);
  Write ('Y0 = [left_y, right_y]                       left_y = ');
  Readln (w);
  Str (w:26, s);
  Y0.a:=left_read(s);
  Write ('                                            right_y = ');
  Readln (w);
  Str (w:26, s);
  Y0.b:=right_read(s);
  if (method<>1) and (method<>9) and (method<>15)
    then begin
           Write ('Y1 = [left_y, right_y]                       left_y = ');
           Readln (w);
           Str (w:26, s);
           Y1.a:=left_read(s);
           Write ('                                            right_y = ');
           Readln (w);
           Str (w:26, s);
           Y1.b:=right_read(s)
         end;
  if (method>=5) and (method<=8) or (method=3) or (method>=19)
    then begin
           Write ('Y2 = [left_y, right_y]                       left_y = ');
           Readln (w);
           Str (w:26, s);
           Y2.a:=left_read(s);
           Write ('                                            right_y = ');
           Readln (w);
           Str (w:26, s);
           Y2.b:=right_read(s)
         end;
  if (method=7) or (method=8)
    then begin
           Write ('Y3 = [left_y, right_y]                       left_y = ');
           Readln (w);
           Str (w:26, s);
           Y3.a:=left_read(s);
           Write ('                                            right_y = ');
           Readln (w);
           Str (w:26, s);
           Y3.b:=right_read(s)
         end;
  find:=False;
  repeat
    Write ('maximum nunmber of steps (m*h <= a)               m = ');
    Readln (m);
    if (m*real_h<=amh) and (m>0)
      then find:=True
  until find;
  repeat
    Write ('output every k steps (k <= m)           output_step = ');
    Readln (output_step)
  until output_step<=m;
  Writeln;
  repeat
    Write ('Output (s - screen, f - file) : ');
    Readln (z)
  until (z='s') or (z='f');
  if z='f'
    then begin
           case method of
              1 : strg:='AB1';
              2 : strg:='N1';
              3 : strg:='AB2';
              4 : strg:='N2';
              5 : strg:='AB3';
              6 : strg:='N3';
              7 : strg:='AB4';
              8 : strg:='N4';
              9 : strg:='AM1';
             10 : strg:='MS1';
             11 : strg:='AM2';
             12 : strg:='MS2';
             13 : strg:='AM3';
             14 : strg:='MS3';
             15 : strg:='ABM1';
             16 : strg:='NMS1';
             17 : strg:='ABM2';
             18 : strg:='NMS2';
             19 : strg:='ABM3';
             20 : strg:='NMS3'
           end;
           strg:='IMM'+strg+'RESULTS.TXT';
           AssignFile (plik, strg);
           Rewrite (plik);
           Writeln ('Output to '+strg+' in current directory')
         end
    else Writeln;
  case method of
     1 : strg:='k = 1, Adams-Bashforth (explicit)';
     2 : strg:='k = 1, Nystrom (explicit)';
     3 : strg:='k = 2, Adams-Bashforth (explicit)';
     4 : strg:='k = 2, Nystrom (explicit, in conventional case the same as for '
               +'k=1)';
     5 : strg:='k = 3, Adams-Bashforth (explicit)';
     6 : strg:='k = 3, Nystrom (explicit)';
     7 : strg:='k = 4, Adams-Bashforth (explicit)';
     8 : strg:='k = 4, Nystrom (explicit)';
     9 : strg:='k = 1, Adams-Moulton (implicit)';
    10 : strg:='k = 1, Milne-Simpson (implicit)';
    11 : strg:='k = 2, Adams-Moulton (implicit)';
    12 : strg:='k = 2, Milne-Simpson (implicit)';
    13 : strg:='k = 3, Adams-Moulton (implicit)';
    14 : strg:='k = 3, Milne-Simpson (implicit, in conventional '
               +'case the same as for k=2)';
    15 : strg:='k = 1, Adams-Bashforth-Moulton predictor-corrector';
    16 : strg:='k = 1, Nystrom-Milne-Simpson predictor-corrector';
    17 : strg:='k = 2, Adams-Bashforth-Moulton predictor-corrector';
    18 : strg:='k = 2, Nystrom-Milne-Simpson predictor-corrector';
    19 : strg:='k = 3, Adams-Bashforth-Moulton predictor-corrector';
    20 : strg:='k = 3, Nystrom-Milne-Simpson predictor-corrector'
  end;
  if z='f'
    then begin
           Writeln (plik, '   ');
           Writeln (plik, strg)
         end
    else Writeln (strg);
  strg:=' accuracy not achieved after 20 iterations';
  case method of
     1         : for n:=1 to m do
                   begin
                     AB1 (T0, Y0, Fdelta, real_h, FTY, PSI1, T1, Y1);
                     if (n mod output_step)=0
                       then Results (T1, Y1, 'e');
                     T0:=T1;
                     Y0:=Y1
                  end;
     2, 3, 4   : for n:=2 to m do
                   begin
                     if method=2
                       then N1 (T0, Y0, Y1, Fdelta, real_h, FTY, PSI1,
                                T1, T2, Y2)
                       else if method=3
                              then AB2 (T0, Y0, Y1, Fdelta, real_h, FTY, PSI1,
                                        T1, T2, Y2)
                              else N2 (T0, Y0, Y1, Fdelta, real_h, FTY, PSI1,
                                       T1, T2, Y2);
                     if (n mod output_step)=0
                       then Results (T2, Y2, 'e');
                     T0:=T1;
                     Y0:=Y1;
                     T1:=T2;
                     Y1:=Y2
                   end;
     5, 6      : for n:=3 to m do
                   begin
                     if method=5
                       then AB3 (T0, Y0, Y1, Y2, Fdelta, real_h, FTY, PSI1,
                                 T1, T2, T3, Y3)
                       else N3 (T0, Y0, Y1, Y2, Fdelta, real_h, FTY, PSI1,
                                T1, T2, T3, Y3);
                     if (n mod output_step)=0
                       then Results (T3, Y3, 'e');
                     T0:=T1;
                     Y0:=Y1;
                     T1:=T2;
                     Y1:=Y2;
                     T2:=T3;
                     Y2:=Y3
                   end;
     7, 8      : for n:=4 to m do
                   begin
                     if method=7
                       then AB4 (T0, Y0, Y1, Y2, Y3, Fdelta, real_h, FTY, PSI1,
                                 T1, T2, T3, T4, Y4)
                       else N4 (T0, Y0, Y1, Y2, Y3, Fdelta, real_h, FTY, PSI1,
                                T1, T2, T3, T4, Y4);
                     if (n mod output_step)=0
                       then Results (T4, Y4, 'e');
                     T0:=T1;
                     Y0:=Y1;
                     T1:=T2;
                     Y1:=Y2;
                     T2:=T3;
                     Y2:=Y3;
                     T3:=T4;
                     Y3:=Y4
                   end;
     9         : begin
                   n:=0;
                   repeat
                     n:=n+1;
                     AM1 (T0, Y0, Fdelta, real_h, FTY, PSI2, find, T1, Y1, l);
                     if not find
                       then Writeln ('For n = ', n, strg)
                       else if (n mod output_step)=0
                              then Results (T1, Y1, 'i');
                     T0:=T1;
                     Y0:=Y1
                   until (n=m) or not find
                 end;
    10, 11, 12 : begin
                   n:=1;
                   repeat
                     n:=n+1;
                     if method=10
                       then MS1 (T0, Y0, Y1, Fdelta, real_h, FTY, PSI2, find,
                                 T1, T2, Y2, l)
                       else if method=11
                              then AM2 (T0, Y0, Y1, Fdelta, real_h, FTY, PSI2,
                                        find, T1, T2, Y2, l)
                              else MS2 (T0, Y0, Y1, Fdelta, real_h, FTY, PSI2,
                                        find, T1, T2, Y2, l);
                     if not find
                       then Writeln ('For n = ', n, strg)
                       else if (n mod output_step)=0
                              then Results (T2, Y2, 'i');
                     T0:=T1;
                     Y0:=Y1;
                     T1:=T2;
                     Y1:=Y2
                   until (n=m) or not find
                 end;
    13         : begin
                   n:=2;
                   repeat
                     n:=n+1;
                     AM3 (T0, Y0, Y1, Y2, Fdelta, real_h, FTY, PSI2, find,
                          T1, T2, T3, Y3, l);
                     if not find
                       then Writeln ('For n = ', n, strg)
                       else if (n mod output_step)=0
                              then Results (T3, Y3, 'i');
                     T0:=T1;
                     Y0:=Y1;
                     T1:=T2;
                     Y1:=Y2;
                     T2:=T3;
                     Y2:=Y3
                   until (n=m) or not find
                 end;
    14         : begin
                   n:=1;
                   repeat
                     n:=n+1;
                     MS3 (T0, Y0, Y1, Fdelta, real_h, FTY, PSI2, find,
                          T1, T2, Y2, l);
                     if not find
                       then Writeln ('For n = ', n, strg)
                       else if (n mod output_step)=0
                              then Results (T2, Y2, 'i');
                     T0:=T1;
                     Y0:=Y1;
                     T1:=T2;
                     Y1:=Y2
                   until (n=m) or not find
                 end;
    15         : begin
                   n:=0;
                   repeat
                     n:=n+1;
                     ABM1 (T0, Y0, Fdelta, real_h, FTY, PSI1, PSI2, find,
                           T1, Y1, l);
                     if not find
                       then Writeln ('For n = ', n, strg)
                       else if (n mod output_step)=0
                              then Results (T1, Y1, 'i');
                     T0:=T1;
                     Y0:=Y1
                   until (n=m) or not find
                 end;
    16, 17, 18 : begin
                   n:=1;
                   repeat
                     n:=n+1;
                     if method=16
                       then NMS1 (T0, Y0, Y1, Fdelta, real_h, FTY, PSI1, PSI2,
                                  find, T1, T2, Y2, l)
                       else if method=17
                              then ABM2 (T0, Y0, Y1, Fdelta, real_h, FTY, PSI1,
                                         PSI2, find, T1, T2, Y2, l)
                              else NMS2 (T0, Y0, Y1, Fdelta, real_h, FTY, PSI1,
                                         PSI2, find, T1, T2, Y2, l);
                     if not find
                       then Writeln ('For n = ', n, strg)
                       else if (n mod output_step)=0
                              then Results (T2, Y2, 'i');
                     T0:=T1;
                     Y0:=Y1;
                     T1:=T2;
                     Y1:=Y2
                   until (n=m) or not find
                 end;
    19, 20     : begin
                   n:=2;
                   repeat
                     n:=n+1;
                     if method=19
                       then ABM3 (T0, Y0, Y1, Y2, Fdelta, real_h, FTY, PSI1,
                                  PSI2, find, T1, T2, T3, Y3, l)
                       else  NMS3 (T0, Y0, Y1, Y2, Fdelta, real_h, FTY, PSI1,
                                   PSI2, find, T1, T2, T3, Y3, l);
                     if not find
                       then Writeln ('For n = ', n, strg)
                       else if (n mod output_step)=0
                              then Results (T3, Y3, 'i');
                     T0:=T1;
                     Y0:=Y1;
                     T1:=T2;
                     Y1:=Y2;
                     T2:=T3;
                     Y2:=Y3
                   until (n=m) or not find
                 end
  end;
  FreeLibrary (ident);
  if z='f'
    then CloseFile (plik);
  Writeln;
  Writeln ('The end');
  Readln
end.