program IPCM_Example_1_2_3;
{$APPTYPE CONSOLE}
uses System.SysUtils, IntervalArithmetic32and64;

// The program solves the initial value problem y' = 0.5y, y(0) = 1 by
// one-, two-, tree- and four-step explicit interval methods of Adams-Bashforth
// type and Nystrom type, by one-, two-, and tree-step implicit interval methods
// of Adams-Moulton type and Milne-Simpson type, and by one-, two, and tree-step
// interval predictor-corrector methods of Adams-Bashforth-Moulton type
// and Nystrom-Milne-Simpson type. It is assumed h = 0.0005, and all starting
// intervals are defined within the program.
// You can direct output to the screen (s) or to the file (f).
// Directing output to the screen, press Enter after each result displayed.

// Important note: Delphi's 64 bit compiler do not support 80-bit Extended
// floating point values on Win64 (Extended = Double on Win64).
// The uTExtendedX87 unit provides for Win64 a replacement FPU-backed 80-bit
// Extended floating point type called TExtendedX87 (see the comments in
// IntervalArithmetic32and64 unit for details).

// Note: Do not compile this program in any Delphi environment older than XE!

function Psi1Explicit (const T, Y : interval) : interval;
begin
  Result:=Y/4
end {Psi1Explicit};

function Psi2Explicit (const T, Y : interval) : interval;
begin
  Result:=Y/8
end {Psi2Explicit};

function Psi3Explicit (const T, Y : interval) : interval;
begin
  Result:=Y/16
end {Psi3Explicit};

function Psi4Explicit (const T, Y : interval) : interval;
begin
  Result:=Y/32
end {Psi4Explicit};

function Psi1Implicit (const T, Y : interval) : interval;
begin
  Result:=Psi2Explicit(T,Y)
end {Psi1Implicit};

function Psi2Implicit (const T, Y : interval) : interval;
begin
  Result:=Psi3Explicit(T,Y)
end {Psi2Implicit};

function Psi3Implicit (const T, Y : interval) : interval;
begin
  Result:=Psi4Explicit(T,Y)
end {Psi3Implicit};

function FTY (const T, Y : interval) : interval;
begin
  Result:=Y/2
end {FTY};

var l, n, m                            : Integer;
    w                                  : Extended;
    delta_t, delta_y, i1, i2, i3, i4,
    i5, h, hh, hh1, Fdelta, T, Y, T0,
    Y0, T1, Y1, T2, Y2, T3, Y3, T4, Y4 : Interval;
    left, right, s                     : string;
    find                               : Boolean;
    kind, z                            : Char;
    plik                               : TextFile;

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
  repeat
    Write ('Output (s - screen, f - file) : ');
    Readln (z)
  until (z='s') or (z='f');
  if z='f'
    then begin
           AssignFile (plik, 'IPCM123RESULTS.TXT');
           Rewrite (plik);
           Writeln ('Output to IPCM123RESULTS.TXT in current directory')
         end
    else Writeln;
  m:=2000;
  s:='0.0005';
  h.a:=left_read(s);
  h.b:=right_read(s);
  delta_t.a:=0;
  delta_t.b:=1;
  delta_y.a:=1;
  s:='1.65';
  delta_y.b:=right_read(s);
  Fdelta:=FTY(delta_t,delta_y);
  if z='s'
    then begin
           Writeln ('Initial value problem: y'' = 0.5y, y(0) = 1');
           Writeln;
           Writeln ('Exact solution:');
           Writeln ('y(0.2) = ', Exp(0.1):26);
           Writeln ('y(0.6) = ', Exp(0.3):26);
           Writeln ('y(1.0) = ', Exp(0.5):26);
           Writeln;
           Readln
         end
    else begin
           Writeln (plik, 'Initial value problem: y'' = 0.5y, y(0) = 1');
           Writeln (plik, ' ');
           Writeln (plik, 'Exact solution:');
           Writeln (plik, '===============');
           Writeln (plik, 'y(0.2) = ', Exp(0.1):26);
           Writeln (plik, 'y(0.6) = ', Exp(0.3):26);
           Writeln (plik, 'y(1.0) = ', Exp(0.5):26)
         end;

//  k=1, Adams-Bashforth (explicit)
  T0:=0;
  Y0:=1;
  s:='0.0005';
  hh.a:=0;
  hh.b:=right_read(s);
  if z='f'
    then begin
           Writeln (plik, '   ');
           Writeln (plik, 'k = 1, Adams-Bashforth (explicit)');
           Writeln (plik, '=================================')
         end
    else Writeln ('k = 1, Adams-Bashforth (explicit)');
  for n:=1 to m do
    begin
      T1:=T0+h;
      T:=T0+hh;
      Y:=Y0+hh*Fdelta;
      i1:=Psi1Explicit(T,Y);
      i1:=h*h*i1;
      i2:=FTY(T0,Y0);
      i2:=h*i2;
      i2:=i2+i1/2;
      Y1:=Y0+i2;
      if (n=400) or (n=1200) or (n=2000)
        then Results (T1, Y1, 'e');
      T0:=T1;
      Y0:=Y1
    end;
  if z='s'
    then Readln;

//  k=1, Nystrom (explicit)
  T0:=0;
  Y0:=1;
  T1:=T0+h;
  s:='1.0002500312526043';
  Y1.a:=left_read(s);
  s:='1.0002500312526044';
  Y1.b:=right_read(s);
  s:='0.0005';
  hh.a:=0;
  hh.b:=right_read(s);
  if z='f'
    then begin
           Writeln (plik, '   ');
           Writeln (plik, 'k = 1, Nystrom (explicit)');
           Writeln (plik, '=========================')
         end
    else Writeln ('k = 1, Nystrom (explicit)');
  for n:=2 to m do
    begin
      T2:=T1+h;
      T:=T1+hh;
      Y:=Y1+hh*Fdelta;
      i1:=Psi1Explicit(T,Y);
      i1:=h*h*(i1-i1);
      i2:=FTY(T1,Y1);
      i2:=2*h*i2;
      i2:=i2+i1/2;
      Y2:=Y0+i2;
      if (n=400) or (n=1200) or (n=2000)
        then Results (T2, Y2, 'e');
      T0:=T1;
      Y0:=Y1;
      T1:=T2;
      Y1:=Y2
    end;
  if z='s'
    then Readln;

//  k=2, Adams-Bashforth (explicit)
  T0:=0;
  Y0:=1;
  T1:=T0+h;
  s:='1.0002500312526043';
  Y1.a:=left_read(s);
  s:='1.0002500312526044';
  Y1.b:=right_read(s);
  s:='-0.0005';
  hh.a:=left_read(s);
  s:='0.0005';
  hh.b:=right_read(s);
  if z='f'
    then begin
           Writeln (plik, '   ');
           Writeln (plik, 'k = 2, Adams-Bashforth (explicit)');
           Writeln (plik, '=================================')
         end
    else Writeln ('k = 2, Adams-Bashforth (explicit)');
  for n:=2 to m do
    begin
      T2:=T1+h;
      T:=T1+hh;
      Y:=Y1+hh*Fdelta;
      i1:=Psi2Explicit(T,Y);
      i2:=h*h*h;
      i2:=5*i2*i1;
      i1:=i2/12;
      i2:=FTY(T1,Y1);
      i3:=FTY(T0,Y0);
      i2:=3*i2-i3;
      i2:=h*i2/2;
      Y2:=Y1+(i2+i1);
      if (n=400) or (n=1200) or (n=2000)
        then Results (T2, Y2, 'e');
      T0:=T1;
      Y0:=Y1;
      T1:=T2;
      Y1:=Y2
    end;
  if z='s'
    then Readln;

//  k=2, Nystrom (explicit, in conventional case the same as for k=1)
  T0:=0;
  Y0:=1;
  T1:=iadd(T0,h);
  s:='1.0002500312526043';
  Y1.a:=left_read(s);
  s:='1.0002500312526044';
  Y1.b:=right_read(s);
  s:='-0.0005';
  hh.a:=left_read(s);
  s:='0.0005';
  hh.b:=right_read(s);
  if z='f'
    then begin
           Writeln (plik, '   ');
           Writeln (plik, 'k = 2, Nystrom (explicit, in conventional case the '
                    +'same as for k=1)');
           Writeln (plik, '==================================================='
                    +'================')
         end
    else Writeln ('k = 2, Nystrom (explicit, in conventional case the same as '
                  +'for k=1)');
  for n:=2 to m do
    begin
      T2:=T1+h;
      T:=T1+hh;
      Y:=Y1+hh*Fdelta;
      i1:=Psi2Explicit(T,Y);
      i1:=5*i1-i1;
      i2:=h*h*h;
      i2:=i2*i1/12;
      i1:=FTY(T1,Y1);
      i1:=2*h*i1;
      Y2:=Y0+(i1+i2);
      if (n=400) or (n=1200) or (n=2000)
        then Results (T2, Y2, 'e');
      T0:=T1;
      Y0:=Y1;
      T1:=T2;
      Y1:=Y2
    end;
  if z='s'
    then Readln;

//  k=3, Adams-Bashforth (explicit)
  T0:=0;
  Y0:=1;
  T1:=T0+h;
  s:='1.0002500312526043';
  Y1.a:=left_read(s);
  s:='1.0002500312526044';
  Y1.b:=right_read(s);
  T2:=T1+h;
  s:='1.0005001250208359';
  Y2.a:=left_read(s);
  s:='1.0005001250208360';
  Y2.b:=right_read(s);
  s:='-0.001';
  hh.a:=left_read(s);
  s:='0.0005';
  hh.b:=right_read(s);
  if z='f'
    then begin
           Writeln (plik, '   ');
           Writeln (plik, 'k = 3, Adams-Bashforth (explicit)');
           Writeln (plik, '=================================')
         end
    else Writeln ('k = 3, Adams-Bashforth (explicit)');
  for n:=3 to m do
    begin
      T3:=T2+h;
      T:=T2+hh;
      Y:=Y2+hh*Fdelta;
      i1:=Psi3Explicit(T,Y);
      i2:=(h*h)*(h*h);
      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:=h*i2/12;
      Y3:=Y2+(i2+i1);
      if (n=400) or (n=1200) or (n=2000)
        then Results (T3, Y3, 'e');
      T0:=T1;
      Y0:=Y1;
      T1:=T2;
      Y1:=Y2;
      T2:=T3;
      Y2:=Y3
    end;
  if z='s'
    then Readln;

//  k=3, Nystrom (explicit)
  T0:=0;
  Y0:=1;
  T1:=T0+h;
  s:='1.0002500312526043';
  Y1.a:=left_read(s);
  s:='1.0002500312526044';
  Y1.b:=right_read(s);
  T2:=T1+h;
  s:='1.0005001250208359';
  Y2.a:=left_read(s);
  s:='1.0005001250208360';
  Y2.b:=right_read(s);
  s:='-0.001';
  hh.a:=left_read(s);
  s:='0.0005';
  hh.b:=right_read(s);
  if z='f'
    then begin
           Writeln (plik, '   ');
           Writeln (plik, 'k = 3, Nystrom (explicit)');
           Writeln (plik, '=========================')
         end
    else Writeln ('k = 3, Nystrom (explicit)');
  for n:=3 to m do
    begin
      T3:=T2+h;
      T:=T2+hh;
      Y:=Y2+hh*Fdelta;
      i1:=Psi3Explicit(T,Y);
      i1:=9*i1-i1;
      i2:=(h*h)*(h*h);
      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:=h*i2/3;
      Y3:=Y1+(i2+i1);
      if (n=400) or (n=1200) or (n=2000)
        then Results (T3, Y3, 'e');
      T0:=T1;
      Y0:=Y1;
      T1:=T2;
      Y1:=Y2;
      T2:=T3;
      Y2:=Y3
    end;
  if z='s'
    then Readln;

//  k=4, Adams-Bashforth (explicit)
  T0:=0;
  Y0:=1;
  T1:=T0+h;
  s:='1.0002500312526043';
  Y1.a:=left_read(s);
  s:='1.0002500312526044';
  Y1.b:=right_read(s);
  T2:=T1+h;
  s:='1.0005001250208359';
  Y2.a:=left_read(s);
  s:='1.0005001250208360';
  Y2.b:=right_read(s);
  T3:=T2+h;
  s:='1.0007502813203256';
  Y3.a:=left_read(s);
  s:='1.0007502813203257';
  Y3.b:=right_read(s);
  s:='-0.0015';
  hh.a:=left_read(s);
  s:='0.0005';
  hh.b:=right_read(s);
  if z='f'
    then begin
           Writeln (plik, '   ');
           Writeln (plik, 'k = 4, Adams-Bashforth (explicit)');
           Writeln (plik, '=================================')
         end
    else Writeln ('k = 4, Adams-Bashforth (explicit)');
  for n:=4 to m do
    begin
      T4:=T3+h;
      T:=T3+hh;
      Y:=Y3+hh*Fdelta;
      i1:=Psi4Explicit(T,Y);
      i2:=(h*h)*(h*h);
      i2:=i2*h;
      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:=h*i2/24;
      Y4:=Y3+(i2+i1);
      if (n=400) or (n=1200) or (n=2000)
        then Results (T4, Y4, 'e');
      T0:=T1;
      Y0:=Y1;
      T1:=T2;
      Y1:=Y2;
      T2:=T3;
      Y2:=Y3;
      T3:=T4;
      Y3:=Y4
    end;
  if z='s'
    then Readln;

//  k=4, Nystrom (explicit)
  T0:=0;
  Y0:=1;
  T1:=T0+h;
  s:='1.0002500312526043';
  Y1.a:=left_read(s);
  s:='1.0002500312526044';
  Y1.b:=right_read(s);
  T2:=T1+h;
  s:='1.0005001250208359';
  Y2.a:=left_read(s);
  s:='1.0005001250208360';
  Y2.b:=right_read(s);
  T3:=T2+h;
  s:='1.0007502813203256';
  Y3.a:=left_read(s);
  s:='1.0007502813203257';
  Y3.b:=right_read(s);
  s:='-0.0015';
  hh.a:=left_read(s);
  s:='0.0005';
  hh.b:=right_read(s);
  if z='f'
    then begin
           Writeln (plik, '   ');
           Writeln (plik, 'k = 4, Nystrom (explicit)');
           Writeln (plik, '=========================')
         end
    else Writeln ('k = 4, Nystrom (explicit)');
  for n:=4 to m do
    begin
      T4:=T3+h;
      T:=T3+hh;
      Y:=Y3+hh*Fdelta;
      i1:=Psi4Explicit(T,Y);
      i1:=251*i1-19*i1;
      i2:=(h*h)*(h*h);
      i2:=i2*h*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:=h*i2/3;
      Y4:=Y2+(i2+i1);
      if (n=400) or (n=1200) or (n=2000)
        then Results (T4, Y4, 'e');
      T0:=T1;
      Y0:=Y1;
      T1:=T2;
      Y1:=Y2;
      T2:=T3;
      Y2:=Y3;
      T3:=T4;
      Y3:=Y4
    end;
  if z='s'
    then Readln;

// k=1, Adams-Moulton (implicit)
  T0:=0;
  Y0:=1;
  s:='-0.0005';
  hh.a:=left_read(s);
  hh.b:=0;
  if z='f'
    then begin
           Writeln (plik, '   ');
           Writeln (plik, 'Adams-Moulton (implicit)');
           Writeln (plik, '========================')
         end
    else Writeln ('k = 1, Adams-Moulton (implicit)');
  n:=0;
  repeat
    n:=n+1;
    T1:=T0+h;
    Y1:=Y0;
    T:=T1+hh;
//    iterations
    l:=0;  // number of iterations
    repeat
      l:=l+1;
      find:=True;
      Y:=Y1+hh*Fdelta;
      i1:=Psi1Implicit(T,Y);
      i2:=h*h*h;
      i1:=i2*i1/12;
      i2:=FTY(T1,Y1);
      i3:=FTY(T0,Y0);
      i2:=i2+i3;
      i2:=h*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
    if not find
      then Writeln ('For n = ', n, ' accuracy not achieved after 20 '
                    +'iterations')
      else if (n=400) or (n=1200) or (n=2000)
             then Results (T1, Y1, 'i');
    T0:=T1;
    Y0:=Y1
  until (n=m) or not find;
  if z='s'
    then Readln;

// k=1, Milne-Simpson (implicit)
  T0:=0;
  Y0:=1;
  T1:=T0+h;
  s:='1.0002500312526043';
  Y1.a:=left_read(s);
  s:='1.0002500312526044';
  Y1.b:=right_read(s);
  s:='-0.0005';
  hh.a:=left_read(s);
  hh.b:=0;
  if z='f'
    then begin
           Writeln (plik, '   ');
           Writeln (plik, 'k = 1, Milne-Simpson (implicit)');
           Writeln (plik, '===============================')
         end
    else Writeln ('k = 1, Milne-Simpson (implicit)');
  n:=1;
  repeat
    n:=n+1;
    T2:=T1+h;
    Y2:=Y1;
    T:=T2+hh;
//    iterations
    l:=0;  // number of iterations
    repeat
      l:=l+1;
      find:=True;
      Y:=Y2+hh*Fdelta;
      i1:=Psi1Implicit(T,Y);
      i1:=5*i1-i1;
      i2:=h*h*h;
      i1:=i2*i1/12;
      i2:=FTY(T1,Y1);
      i2:=2*h*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
    if not find
      then Writeln ('For n = ', n, ' accuracy not achieved after 20 '
                    +'iterations')
      else if (n=400) or (n=1200) or (n=2000)
             then Results (T2, Y2, 'i');
    T0:=T1;
    Y0:=Y1;
    T1:=T2;
    Y1:=Y2
  until (n=m) or not find;
  if z='s'
    then Readln;

// k=2, Adams-Moulton (implicit)
  T0:=0;
  Y0:=1;
  T1:=T0+h;
  s:='1.0002500312526043';
  Y1.a:=left_read(s);
  s:='1.0002500312526044';
  Y1.b:=right_read(s);
  s:='-0.001';
  hh.a:=left_read(s);
  hh.b:=0;
  if z='f'
    then begin
           Writeln (plik, '   ');
           Writeln (plik, 'k = 2, Adams-Moulton (implicit)');
           Writeln (plik, '===============================')
         end
    else Writeln ('k = 2, Adams-Moulton (implicit)');
  n:=1;
  repeat
    n:=n+1;
    T2:=T1+h;
    Y2:=Y1;
    T:=T2+hh;
//    iterations
    l:=0;  // number of iterations
    repeat
      l:=l+1;
      find:=True;
      Y:=Y2+hh*Fdelta;
      i1:=Psi2Implicit(T,Y);
      i2:=(h*h)*(h*h);
      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:=h*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
    if not find
      then Writeln ('For n = ', n, ' accuracy not achieved after 20 '
                    +'iterations')
      else if (n=400) or (n=1200) or (n=2000)
             then Results (T2, Y2, 'i');
    T0:=T1;
    Y0:=Y1;
    T1:=T2;
    Y1:=Y2
  until (n=m) or not find;
  if z='s'
    then Readln;

// k=2, Milne-Simpson (implicit)
  T0:=0;
  Y0:=1;
  T1:=T0+h;
  s:='1.0002500312526043';
  Y1.a:=left_read(s);
  s:='1.0002500312526044';
  Y1.b:=right_read(s);
  s:='-0.001';
  hh.a:=left_read(s);
  hh.b:=0;
  if z='f'
    then begin
           Writeln (plik, '   ');
           Writeln (plik, 'k = 2, Milne-Simpson (implicit)');
           Writeln (plik, '===============================')
         end
    else Writeln ('k = 2, Milne-Simpson (implicit)');
  n:=1;
  repeat
    n:=n+1;
    T2:=T1+h;
    Y2:=Y1;
    T:=T2+hh;
//    iterations
    l:=0;  // number of iterations
    repeat
      l:=l+1;
      find:=True;
      Y:=Y2+hh*Fdelta;
      i1:=Psi2Implicit(T,Y);
      i1:=i1-i1;
      i2:=(h*h)*(h*h);
      i1:=i2*i1/24;
      i2:=FTY(T2,Y2);
      i3:=FTY(T1,Y1);
      i2:=i2+4*i3;
      i3:=FTY(T0,Y0);
      i2:=i2+i3;
      i2:=h*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
    if not find
      then Writeln ('For n = ', n, ' accuracy not achieved after 20 '
                    +'iterations')
      else if (n=400) or (n=1200) or (n=2000)
             then Results (T2, Y2, 'i');
    T0:=T1;
    Y0:=Y1;
    T1:=T2;
    Y1:=Y2
  until (n=m) or not find;
  if z='s'
    then Readln;

// k=3, Adams-Moulton (implicit)
  T0:=0;
  Y0:=1;
  T1:=T0+h;
  s:='1.0002500312526043';
  Y1.a:=left_read(s);
  s:='1.0002500312526044';
  Y1.b:=right_read(s);
  T2:=T1+h;
  s:='1.0005001250208359';
  Y2.a:=left_read(s);
  s:='1.0005001250208360';
  Y2.b:=right_read(s);
  s:='-0.0015';
  hh.a:=left_read(s);
  hh.b:=0;
  if z='f'
    then begin
           Writeln (plik, '   ');
           Writeln (plik, 'k = 3, Adams-Moulton (implicit)');
           Writeln (plik, '===============================')
         end
    else Writeln ('k = 3, Adams-Moulton (implicit)');
  n:=2;
  repeat
    n:=n+1;
    T3:=T2+h;
    Y3:=Y2;
    T:=T3+hh;
//    iterations
    l:=0;  // number of iterations
    repeat
      l:=l+1;
      find:=True;
      Y:=Y3+hh*Fdelta;
      i1:=Psi3Implicit(T,Y);
      i2:=h*h*h;
      i2:=i2*(h*h);
      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:=h*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
    if not find
      then Writeln ('For n = ', n, ' accuracy not achieved after 20 '
                    +'iterations')
      else if (n=400) or (n=1200) or (n=2000)
             then Results (T3, Y3, 'i');
    T0:=T1;
    Y0:=Y1;
    T1:=T2;
    Y1:=Y2;
    T2:=T3;
    Y2:=Y3
  until (n=m) or not find;
  if z='s'
    then Readln;

// k=3, Milne-Simpson (implicit, in conventional case the same as for k=2)
  T0:=0;
  Y0:=1;
  T1:=T0+h;
  s:='1.0002500312526043';
  Y1.a:=left_read(s);
  s:='1.0002500312526044';
  Y1.b:=right_read(s);
  s:='-0.0015';
  hh.a:=left_read(s);
  hh.b:=0;
  if z='f'
    then begin
           Writeln (plik, '   ');
           Writeln (plik, 'k = 3, Milne-Simpson (implicit, in conventional '
                    +'case the same as for k=2)');
           Writeln (plik, '================================================'
                    +'=========================')
         end
    else begin
           Writeln ('k = 3, Milne-Simpson (implicit, in conventional case '
                    +'the same as for k=2)')
         end;
  n:=1;
  repeat
    n:=n+1;
    T2:=T1+h;
    Y2:=Y1;
    T:=T2+hh;
//    iterations
    l:=0;  // number of iterations
    repeat
      l:=l+1;
      find:=True;
      Y:=Y2+hh*Fdelta;
      i1:=Psi3Implicit(T,Y);
      i1:=11*i1-19*i1;
      i2:=(h*h)*(h*h);
      i2:=i2*h;
      i1:=i2*i1/720;
      i2:=FTY(T2,Y2);
      i3:=FTY(T1,Y1);
      i2:=i2+4*i3;
      i3:=FTY(T0,Y0);
      i2:=i2+i3;
      i2:=h*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
    if not find
      then Writeln ('For n = ', n, ' accuracy not achieved after 20 '
                    +'iterations')
      else if (n=400) or (n=1200) or (n=2000)
             then Results (T2, Y2, 'i');
    T0:=T1;
    Y0:=Y1;
    T1:=T2;
    Y1:=Y2
  until (n=m) or not find;
  if z='s'
    then Readln;

// k=1, Adams-Bashforth-Moulton predictor-corrector
  T0:=0;
  Y0:=1;
  s:='0.0005';
  hh.a:=0;
  hh.b:=right_read(s);
  s:='-0.0005';
  hh1.a:=left_read(s);
  hh1.b:=0;
  if z='f'
    then begin
           Writeln (plik, '   ');
           Writeln (plik, 'k = 1, Adams-Bashforth-Moulton predictor-corrector');
           Writeln (plik, '==================================================')
         end
    else Writeln ('k = 1, Adams-Bashforth-Moulton predictor-corrector');
  n:=0;
  repeat
    n:=n+1;
    T1:=T0+h;
    T:=T0+hh;
    Y:=Y0+hh*Fdelta;
    i1:=Psi1Explicit(T,Y);
    i1:=h*h*i1;
    i2:=FTY(T0,Y0);
    i2:=h*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:=Psi1Implicit(T,Y);
      i2:=h*h*h;
      i1:=i2*i1/12;
      i2:=FTY(T1,Y1);
      i3:=FTY(T0,Y0);
      i2:=i2+i3;
      i2:=h*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
    if not find
      then Writeln ('For n = ', n, ' accuracy not achieved after 20 '
                    +'iterations')
      else if (n=400) or (n=1200) or (n=2000)
             then Results (T1, Y1, 'i');
    T0:=T1;
    Y0:=Y1
  until (n=m) or not find;
  if z='s'
    then Readln;

// k=1, Nystrom-Milne-Simpson predictor-corrector
  T0:=0;
  Y0:=1;
  T1:=T0+h;
  s:='1.0002500312526043';
  Y1.a:=left_read(s);
  s:='1.0002500312526044';
  Y1.b:=right_read(s);
  s:='0.0005';
  hh.a:=0;
  hh.b:=right_read(s);
  s:='-0.0005';
  hh1.a:=left_read(s);
  hh1.b:=0;
  if z='f'
    then begin
           Writeln (plik, '   ');
           Writeln (plik, 'k = 1, Nystrom-Milne-Simpson predictor-corrector');
           Writeln (plik, '================================================')
         end
    else Writeln ('k = 1, Nystrom-Milne-Simpson predictor-corrector');
  n:=1;
  repeat
    n:=n+1;
    T2:=T1+h;
    T:=T1+hh;
    Y:=Y1+hh*Fdelta;
    i1:=Psi1Explicit(T,Y);
    i1:=h*h*(i1-i1);
    i2:=FTY(T1,Y1);
    i2:=2*h*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:=Psi1Implicit(T,Y);
      i1:=5*i1-i1;
      i2:=h*h*h;
      i1:=i2*i1/12;
      i2:=FTY(T1,Y1);
      i2:=2*h*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
    if not find
      then Writeln ('For n = ', n, ' accuracy not achieved after 20 '
                    +'iterations')
      else if (n=400) or (n=1200) or (n=2000)
             then Results (T2, Y2, 'i');
    T0:=T1;
    Y0:=Y1;
    T1:=T2;
    Y1:=Y2
  until (n=m) or not find;
  if z='s'
    then Readln;

// k=2, Adams-Bashforth-Moulton predictor-corrector
  T0:=0;
  Y0:=1;
  T1:=T0+h;
  s:='1.0002500312526043';
  Y1.a:=left_read(s);
  s:='1.0002500312526044';
  Y1.b:=right_read(s);
  s:='-0.0005';
  hh.a:=left_read(s);
  s:='0.0005';
  hh.b:=right_read(s);
  s:='-0.001';
  hh1.a:=left_read(s);
  hh1.b:=0;
  if z='f'
    then begin
           Writeln (plik, '   ');
           Writeln (plik, 'k = 2, Adams-Bashforth-Moulton predictor-corrector');
           Writeln (plik, '==================================================')
         end
    else Writeln ('k = 2, Adams-Bashforth-Moulton predictor-corrector');
  n:=1;
  repeat
    n:=n+1;
    T2:=T1+h;
    T:=T1+hh;
    Y:=Y1+hh*Fdelta;
    i1:=Psi2Explicit(T,Y);
    i2:=h*h*h;
    i2:=5*i2*i1;
    i1:=i2/12;
    i2:=FTY(T1,Y1);
    i3:=FTY(T0,Y0);
    i2:=3*i2-i3;
    i2:=h*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:=Psi2Implicit(T,Y);
      i2:=(h*h)*(h*h);
      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:=h*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
    if not find
      then Writeln ('For n = ', n, ' accuracy not achieved after 20 '
                    +'iterations')
      else if (n=400) or (n=1200) or (n=2000)
             then Results (T2, Y2, 'i');
    T0:=T1;
    Y0:=Y1;
    T1:=T2;
    Y1:=Y2
  until (n=m) or not find;
  if z='s'
    then Readln;

// k=2, Nystrom-Milne-Simpson predictor-corrector
  T0:=0;
  Y0:=1;
  T1:=T0+h;
  s:='1.0002500312526043';
  Y1.a:=left_read(s);
  s:='1.0002500312526044';
  Y1.b:=right_read(s);
  s:='-0.0005';
  hh.a:=left_read(s);
  s:='0.0005';
  hh.b:=right_read(s);
  s:='-0.001';
  hh1.a:=left_read(s);
  hh1.b:=0;
  if z='f'
    then begin
           Writeln (plik, '   ');
           Writeln (plik, 'k = 2, Nystrom-Milne-Simpson predictor-corrector');
           Writeln (plik, '================================================')
         end
    else Writeln ('k = 2, Nystrom-Milne-Simpson predictor-corrector');
  n:=1;
  repeat
    n:=n+1;
    T2:=T1+h;
    T:=T1+hh;
    Y:=Y1+hh*Fdelta;
    i1:=Psi2Explicit(T,Y);
    i1:=5*i1-i1;
    i2:=h*h*h;
    i2:=i2*i1/12;
    i1:=FTY(T1,Y1);
    i1:=2*h*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:=Psi2Implicit(T,Y);
      i1:=i1-i1;
      i2:=(h*h)*(h*h);
      i1:=i2*i1/24;
      i2:=FTY(T2,Y2);
      i3:=FTY(T1,Y1);
      i2:=i2+4*i3;
      i3:=FTY(T0,Y0);
      i2:=i2+i3;
      i2:=h*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
    if not find
      then Writeln ('For n = ', n, ' accuracy not achieved after 20 '
                    +'iterations')
      else if (n=400) or (n=1200) or (n=2000)
             then Results (T2, Y2, 'i');
    T0:=T1;
    Y0:=Y1;
    T1:=T2;
    Y1:=Y2
  until (n=m) or not find;
  if z='s'
    then Readln;

// k=3, Adams-Bashforth-Moulton predictor-corrector
  T0:=0;
  Y0:=1;
  T1:=T0+h;
  s:='1.0002500312526043';
  Y1.a:=left_read(s);
  s:='1.0002500312526044';
  Y1.b:=right_read(s);
  T2:=T1+h;
  s:='1.0005001250208359';
  Y2.a:=left_read(s);
  s:='1.0005001250208360';
  Y2.b:=right_read(s);
  s:='-0.001';
  hh.a:=left_read(s);
  s:='0.0005';
  hh.b:=right_read(s);
  s:='-0.0015';
  hh1.a:=left_read(s);
  hh1.b:=0;
  if z='f'
    then begin
           Writeln (plik, '   ');
           Writeln (plik, 'k = 3, Adams-Bashforth-Moulton predictor-corrector');
           Writeln (plik, '==================================================')
         end
    else Writeln ('k = 3, Adams-Bashforth-Moulton predictor-corrector');
  n:=2;
  repeat
    n:=n+1;
    T3:=T2+h;
    T:=T2+hh;
    Y:=Y2+hh*Fdelta;
    i1:=Psi3Explicit(T,Y);
    i2:=(h*h)*(h*h);
    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:=h*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:=Psi3Implicit(T,Y);
      i2:=h*h*h;
      i2:=i2*(h*h);
      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:=h*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
    if not find
      then Writeln ('For n = ', n, ' accuracy not achieved after 20 '
                    +'iterations')
      else if (n=400) or (n=1200) or (n=2000)
             then Results (T3, Y3, 'i');
    T0:=T1;
    Y0:=Y1;
    T1:=T2;
    Y1:=Y2;
    T2:=T3;
    Y2:=Y3
  until (n=m) or not find;
  if z='s'
    then Readln;

// k=3, Nystrom-Milne-simpson predictor-corrector
  T0:=0;
  Y0:=1;
  T1:=T0+h;
  s:='1.0002500312526043';
  Y1.a:=left_read(s);
  s:='1.0002500312526044';
  Y1.b:=right_read(s);
  T2:=T1+h;
  s:='1.0005001250208359';
  Y2.a:=left_read(s);
  s:='1.0005001250208360';
  Y2.b:=right_read(s);
  s:='-0.001';
  hh.a:=left_read(s);
  s:='0.0005';
  hh.b:=right_read(s);
  s:='-0.0015';
  hh1.a:=left_read(s);
  hh1.b:=0;
  if z='f'
    then begin
           Writeln (plik, '   ');
           Writeln (plik, 'k = 3, Nystrom-Milne-Simpson predictor-corrector');
           Writeln (plik, '================================================')
         end
    else Writeln ('k = 3, Nystrom-Milne-Simpson predictor-corrector');
  n:=2;
  repeat
    n:=n+1;
    T3:=T2+h;
    T:=T2+hh;
    Y:=Y2+hh*Fdelta;
    i1:=Psi3Explicit(T,Y);
    i1:=9*i1-i1;
    i2:=(h*h)*(h*h);
    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:=h*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:=Psi3Implicit(T,Y);
      i1:=11*i1-19*i1;
      i2:=(h*h)*(h*h);
      i2:=i2*h;
      i1:=i2*i1/720;
      i2:=FTY(T3,Y3);
      i3:=FTY(T2,Y2);
      i2:=i2+4*i3;
      i3:=FTY(T1,Y1);
      i2:=iadd(i2,i3);
      i2:=h*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
    if not find
      then Writeln ('For n = ', n, ' accuracy not achieved after 20 '
                    +'iterations')
      else if (n=400) or (n=1200) or (n=2000)
             then Results (T3, Y4, 'i');
    T0:=T1;
    Y0:=Y1;
    T1:=T2;
    Y1:=Y2;
    T2:=T3;
    Y2:=Y3
  until (n=m) or not find;
  if z='f'
    then CloseFile (plik);
  Writeln;
  Writeln ('The end. Press Enter ...');
  Readln
end.