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

// The program solves a special case of Hill equations given by the initial
// value problem y'[l] = y[l+2] (l = 1, 2), y'[3] = -(kappa/r^3)y[1],
// y'[4] = -(kappa/r^3)y[2], y[1](0) = y[4](0) = 1, y[2](0) = y[3](0) = 0,
// where r = Sqrt(y[1]^2 + y[2]^2), kappa = 1, by three-step interval
// predictor-corrector method of Nystrom-Milne-Simpson type. It is assumed
// h = 0.0005, and tmin <= t <= tmax, for which you should enter tmax = k,
// where k = 1, ... , 5 (tmin = tmax - 1).
// All starting intervals are defined within the program. The starting intervals
// for k > 1 are taken from the intervals obtained for k - 1.
// 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!

type array4interval = array [1..4] of interval;

function Psi3Explicit (const T : interval;
                       const Y : array4interval) : array4interval;
var psi, r2, r4, r5, r6, r7, r8,
    sqr_y1, sqr_y2, sqr_y3, sqr_y4,
    sqr2_y1, sqr2_y2,
    w1, w2, w3, w4, w5, y1y2, y3y4     : interval;
begin
  sqr_y1:=Y[1]*Y[1];
  sqr_y2:=Y[2]*Y[2];
  r2:=sqr_y1+sqr_y2;
  w1.a:=Sqrt(r2.a);
  w1.b:=Sqrt(r2.b);
  r4:=r2*r2;
  r6:=r4*r2;
  r8:=r4*r4;
  r5:=w1*r4;
  r7:=w1*r6;
  sqr_y3:=Y[3]*Y[3];
  sqr_y4:=Y[4]*Y[4];
  y1y2:=Y[1]*Y[2];
  y3y4:=Y[3]*Y[4];
  w1:=3-5*sqr_y1/r2;
  w2:=1-5*sqr_y1/r2;
  w3:=1-5*sqr_y2/r2;
  w4:=1-3*sqr_y1/r2;
  psi:=(3*Y[1]/r5)*(w1*sqr_y3);
  w1:=(6*Y[2]/r5)*(w2*y3y4);
  psi:=psi+w1;
  w1:=(3*Y[1]/r5)*(w3*sqr_y4);
  psi:=psi+w1;
  w1:=(1/r6)*(w4*Y[1]);
  psi:=psi+w1;
  w1:=(3*y1y2/r8)*Y[2];
  Result[1]:=psi-w1;
  w1:=3-5*sqr_y2/r2;
  w4:=1-3*sqr_y2/r2;
  psi:=(3*Y[2]/r5)*(w1*sqr_y4);
  w1:=(3*Y[2]/r5)*(w2*sqr_y3);
  psi:=psi+w1;
  w1:=(6*Y[1]/r5)*(w3*y3y4);
  psi:=psi+w1;
  w1:=(3*y1y2/r8)*Y[1];
  psi:=psi-w1;
  w1:=(1/r6)*(w4*Y[2]);
  Result[2]:=psi+w1;
  w1:=35*(sqr_y1*sqr_y2)/r4;
  w1:=w1-5*(sqr_y1+sqr_y2)/r2;
  w1:=1+w1;
  w2:=9*(sqr_y1+sqr_y2)/r2;
  w2:=w2-4;
  w3:=15*y1y2/r7;
  w4:=6*y1y2/r8;
  sqr2_y1:=sqr_y1*sqr_y1;
  sqr2_y2:=sqr_y2*sqr_y2;
  w5:=3-30*sqr_y1/r2;
  w5:=w5+35*sqr2_y1/r4;
  psi:=(3/r5)*(w5*(sqr_y3*Y[3]));
  w5:=21*sqr_y1/r2-9;
  w5:=w3*(w5*(sqr_y3*Y[4]));
  psi:=psi+w5;
  w5:=(9/r5)*(w1*(Y[3]*sqr_y4));
  psi:=psi+w5;
  w5:=7*sqr_y2/r2-3;
  w5:=w3*(w5*(sqr_y4*Y[4]));
  psi:=psi+w5;
  w5:=54*(sqr_y1*sqr_y2)/r4;
  w5:=54*sqr2_y1/r4+w5;
  w5:=w5-9*sqr_y2/r2;
  w5:=w5-33*sqr_y1/r2;
  w5:=1+w5;
  w5:=(1/r6)*(w5*Y[3]);
  psi:=psi+w5;
  w5:=w4*(w2*Y[4]);
  Result[3]:=psi+w4;
  w5:=7*sqr_y1/r2-3;
  psi:=w3*w5;
  w5:=(9/r5)*w1;
  w5:=w5*(sqr_y3*Y[4]);
  psi:=psi+w5;
  w5:=21*sqr_y2/r2-9;
  psi:=psi+w3*w5;
  w5:=35*sqr2_y2/r4;
  w5:=w5-30*sqr_y2/r2;
  w5:=(3/r5)*(3+w5);
  psi:=psi+(w5*(sqr_y4*Y[4]));
  w5:=w4*(w2*Y[3]);
  psi:=psi+w5;
  w5:=54*sqr2_y2/r4;
  w5:=w5+54*(sqr_y1*sqr_y2)/r4;
  w5:=w5-33*sqr_y2/r2;
  w5:=w5-9*sqr_y1/r2;
  w5:=((1/r6)*(1+w5))*Y[4];
  Result[4]:=psi+w5
end {Psi3Explicit};

function Psi4Explicit (const T : interval;
                       const Y : array4interval) : array4interval;
var psi, r2, r4, r5, r6, r7, r8, r9,
    sqr_y1, sqr_y2, sqr_y3, sqr_y4,
    sqr2_y1, sqr2_y2, sqr2_y3, sqr2_y4,
    w1, w2, w3, w4, w5, w6, y1y2, y3y4 : interval;
begin
  sqr_y1:=Y[1]*Y[1];
  sqr_y2:=Y[2]*Y[2];
  r2:=sqr_y1+sqr_y2;
  w1.a:=Sqrt(r2.a);
  w1.b:=Sqrt(r2.b);
  r4:=r2*r2;
  r6:=r4*r2;
  r8:=r4*r4;
  r5:=w1*r4;
  r7:=w1*r6;
  r9:=w1*r8;
  sqr_y3:=Y[3]*Y[3];
  sqr_y4:=Y[4]*Y[4];
  y1y2:=Y[1]*Y[2];
  y3y4:=Y[3]*Y[4];
  w1:=35*(sqr_y1*sqr_y2)/r4;
  w1:=w1-5*(sqr_y1+sqr_y2)/r2;
  w1:=1+w1;
  w2:=9*(sqr_y1+sqr_y2)/r2;
  w2:=w2-4;
  w3:=15*y1y2/r7;
  w4:=6*y1y2/r8;
  sqr2_y1:=sqr_y1*sqr_y1;
  sqr2_y2:=sqr_y2*sqr_y2;
  sqr2_y3:=sqr_y3*sqr_y3;
  sqr2_y4:=sqr_y4*sqr_y4;
  w5:=3-30*sqr_y1/r2;
  w5:=w5+35*sqr2_y1/r4;
  psi:=(3/r5)*(w5*(sqr_y3*Y[3]));
  w5:=21*sqr_y1/r2-9;
  w5:=w3*(w5*(sqr_y3*Y[4]));
  psi:=psi+w5;
  w5:=(9/r5)*(w1*(Y[3]*sqr_y4));
  psi:=psi+w5;
  w5:=7*sqr_y2/r2-3;
  w5:=w3*(w5*(sqr_y4*Y[4]));
  psi:=psi+w5;
  w5:=54*(sqr_y1*sqr_y2)/r4;
  w5:=54*sqr2_y1/r4+w5;
  w5:=w5-9*sqr_y2/r2;
  w5:=w5-33*sqr_y1/r2;
  w5:=1+w5;
  w5:=(1/r6)*(w5*Y[3]);
  psi:=psi+w5;
  w5:=w4*(w2*Y[4]);
  Result[1]:=psi+w4;
  w5:=7*sqr_y1/r2-3;
  psi:=w3*(w5*(sqr_y3*Y[3]));
  w5:=(9/r5)*w1;
  w5:=w5*(sqr_y3*Y[4]);
  psi:=psi+w5;
  w5:=21*sqr_y2/r2-9;
  psi:=psi+w3*(w5*(Y[3]*sqr_y4));
  w5:=35*sqr2_y2/r4;
  w5:=w5-30*sqr_y2/r2;
  w5:=(3/r5)*(3+w5);
  psi:=psi+w5*(sqr_y4*Y[4]);
  w5:=w4*(w2*Y[3]);
  psi:=psi+w5;
  w5:=54*sqr2_y2/r4;
  w5:=w5+54*(sqr_y1*sqr_y2)/r4;
  w5:=w5-33*sqr_y2/r2;
  w5:=w5-9*sqr_y1/r2;
  w5:=(1/r6)*(1+w5)*Y[4];
  Result[2]:=psi+w5;
  w1:=14*sqr_y1/r2;
  w1:=w1-21*sqr2_y1/r4;
  w1:=w1-1;
  w2:=14*sqr_y2/r2;
  w2:=w1-21*sqr2_y2/r4;
  w2:=w1-1;
  w3:=7*sqr_y1/r2;
  w3:=w3+21*sqr_y2/r2;
  w4:=63*(sqr_y1*sqr_y2)/r4;
  w3:=w3-w4-3;
  w4:=21*sqr_y1/r2-w4;
  w4:=w4+7*sqr_y2/r2;
  w4:=w4-3;
  w5:=33*(sqr_y1+sqr_y2)/r2;
  w5:=w5-54*(sqr2_y1+sqr2_y2)/r4;
  w5:=w5-108*(sqr_y1*sqr_y2)/r4;
  w5:=w5-1;
  w6:=70*sqr_y1/r2;
  w6:=w6-63*sqr2_y1/r4;
  w6:=(w6-15)*sqr2_y3;
  psi:=(15*Y[1]/r7)*w6;
  w6:=w1*((sqr_y3*Y[3])*Y[4]);
  w6:=(180*Y[2]/r7)*w6;
  psi:=psi+w6;
  w6:=w3*(sqr_y3*sqr_y4);
  w6:=(90*Y[1]/r7)*w6;
  psi:=psi+w6;
  w6:=w4*(Y[3]*(sqr_y4*Y[4]));
  w6:=(60*Y[2]/r7)*w6;
  psi:=psi+w6;
  w6:=w2*sqr2_y4;
  w6:=(45*Y[1]/r7)*w6;
  psi:=psi+w6;
  w6:=(250*sqr_y1+105*sqr_y2)/r2;
  w6:=w6-285*sqr2_y1/r4;
  w6:=w6-285*(sqr_y1*sqr_y2)/r4;
  w6:=(w6-33)*sqr_y3;
  w6:=3*Y[1]*w6/r8;
  psi:=psi+w6;
  w6:=(181*sqr_y1+36*sqr_y2)/r2;
  w6:=w6-285*sqr2_y1/r4;
  w6:=w6-285*(sqr_y1*sqr_y2)/r4;
  w6:=w6*(Y[3]*Y[4]);
  w6:=6*Y[2]*w6/r8;
  psi:=psi+w6;
  w6:=(33*sqr_y1+178*sqr_y2)/r2;
  w6:=w6-285*(sqr_y1*sqr_y2)/r4;
  w6:=w6-285*sqr2_y2/r4;
  w6:=(w6-11)*sqr2_y4;
  w6:=3*Y[1]*w6/r8;
  psi:=psi+w6;
  w6:=Y[1]*w5/r9;
  Result[3]:=psi+w6;
  w6:=w1*sqr2_y3;
  psi:=45*(Y[2]*w6)/r7;
  w6:=w3*((sqr_y3*Y[3])*Y[4]);
  w6:=60*(Y[1]*w6)/r7;
  psi:=psi+w6;
  w6:=w4*(sqr_y3*sqr_y3);
  w6:=90*(Y[2]*w6)/r7;
  psi:=psi+w6;
  w6:=w2*(Y[3]*(sqr_y4*Y[4]));
  w6:=180*(Y[1]*w6)/r7;
  psi:=psi+w6;
  w6:=70*sqr_y2/r2;
  w6:=w6-63*sqr2_y2/r4;
  w6:=(w6-15)*sqr2_y4;
  w6:=15*(Y[2]*w6)/r7;
  psi:=psi+w6;
  w6:=(178*sqr_y1+33*sqr_y2)/r2;
  w6:=w6-285*sqr2_y1/r4;
  w6:=w6-285*(sqr_y1*sqr_y2)/r4;
  w6:=(w6-11)*sqr2_y3;
  w6:=3*(Y[2]*w6)/r8;
  psi:=psi+w6;
  w6:=(36*sqr_y1+181*sqr_y2)/r2;
  w6:=w6-285*(sqr_y1*sqr_y2)/r4;
  w6:=w6-285*sqr2_y2/r4;
  w6:=(w6-11)*(Y[3]*Y[4]);
  w6:=6*(Y[1]*w6)/r8;
  psi:=psi+w6;
  w6:=(105*sqr_y1+250*sqr_y2)/r2;
  w6:=w6-285*(sqr_y1*sqr_y2)/r4;
  w6:=w6-285*sqr2_y2/r4;
  w6:=(w6-33)*sqr2_y4;
  w6:=3*(Y[2]*w6)/r8;
  psi:=psi+w6;
  w6:=Y[2]*w5/r9;
  Result[4]:=psi+w6
end {Psi4Explicit};

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

function FTY (const T : interval;
              const Y : array4interval) : array4interval;
var r, r2 : interval;
begin
  Result[1]:=Y[3];
  Result[2]:=Y[4];
  r2:=Y[1]*Y[1]+Y[2]*Y[2];
  r.a:=Sqrt(r2.a);
  r.b:=Sqrt(r2.b);
  r:=r2*r;
  Result[3]:=-Y[1]/r;
  Result[4]:=-Y[2]/r
end {FTY};

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

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

begin
  Writeln ('The initial value problem y''[l] = y[l+2] (l = 1, 2),');
  Writeln ('                          y''[3] = -y[1]/r^3,');
  Writeln ('                          y''[4] = -y[2]/r^3,');
  Writeln ('                          y[1](0) = y[4](0) = 1,'
           +' y[2](0) = y[3](0) = 0,');
  Writeln ('where r = Sqrt(y[1]^2 + y[2]^2), tmin <= t <= tmax');
  Writeln;
  repeat
    Write ('Enter tmax (tmax must be integer and 1 <= tmax <= 5)   : ');
    Readln (k)
  until (k>=1) and (k<=5);
  repeat
    Write ('Do you need starting intervals for the next tmax (y/n) : ');
    Readln (next)
  until (next='y') or (next='n');
  Writeln;
  repeat
    Write ('Output (s - screen, f - file) : ');
    Readln (z)
  until (z='s') or (z='f');
  if z='f'
    then begin
           AssignFile (plik, 'IPCHILLRESULTS.TXT');
           Rewrite (plik);
           Writeln ('Output to IPCHILLRESULTS.TXT in current directory')
         end
    else Writeln;
  if z='s'
    then begin
           Writeln ('Exact solution:');
           Writeln ('y[1](1.0) = ', Cos(1.0):26,
                    '  y[2](1.0) = ', Sin(1.0):26);
           Writeln ('y[3](1.0) = ', -Sin(1.0):26,
                    '  y[4](1.0) = ', Cos(1.0):26);
           Writeln ('y[1](3.0) = ', Cos(3.0):26,
                    '  y[2](3.0) = ', Sin(3.0):26);
           Writeln ('y[3](3.0) = ', -Sin(3.0):26,
                    '  y[4](3.0) = ', Cos(3.0):26);
           Writeln ('y[1](5.0) = ', Cos(5.0):26,
                    '  y[2](5.0) = ', Sin(5.0):26);
           Writeln ('y[3](5.0) = ', -Sin(5.0):26,
                    '  y[4](5.0) = ', Cos(5.0):26);
           Writeln;
           Readln
         end
    else begin
           Writeln (plik, 'Exact solution:');
           Writeln (plik, '===============');
           Writeln (plik, 'y[1](1.0) = ', Cos(1.0):26,
                    '  y[2](1.0) = ', Sin(1.0):26);
           Writeln (plik, 'y[3](1.0) = ', -Sin(1.0):26,
                    '  y[4](1.0) = ', Cos(1.0):26);
           Writeln (plik, 'y[1](3.0) = ', Cos(3.0):26,
                    '  y[2](3.0) = ', Sin(3.0):26);
           Writeln (plik, 'y[3](3.0) = ', -Sin(3.0):26,
                    '  y[4](3.0) = ', Cos(3.0):26);
           Writeln (plik, 'y[1](5.0) = ', Cos(5.0):26,
                    '  y[2](5.0) = ', Sin(5.0):26);
           Writeln (plik, 'y[3](5.0) = ', -Sin(5.0):26,
                    '  y[4](5.0) = ', Cos(5.0):26)
         end;
  Writeln ('Be patient ... Calculations may take a lot of minutes.');
  Writeln;
  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');
  m:=2002;
  s:='0.0005';
  h.a:=left_read(s);
  h.b:=right_read(s);
  case k of
    1 : begin
          delta_t.a:=0;
          s:='1.001';
        end;
    2 : begin
          delta_t.a:=1;
          s:='2.001'
        end;
    3 : begin
          delta_t.a:=2;
          s:='3.001'
        end;
    4 : begin
          delta_t.a:=3;
          s:='4.001'
        end;
    5 : begin
          delta_t.a:=4;
          s:='5.001'
        end;
  end;
  delta_t.b:=right_read(s);
  if k=4
    then delta_y[1].a:=-1
    else begin
           case k of
             1 : s:='0.4';
             2 : s:='-0.5';
             3 : s:='-1.1';
             5 : s:='-0.6'
           end;
           delta_y[1].a:=left_read(s);
         end;
  case k of
    1 : s:='1.1';
    2 : s:='0.6';
    3 : s:='-0.3';
    4 : s:='-0.6';
    5 : s:='0.4'
  end;
  delta_y[1].b:=right_read(s);
  if k=3
    then delta_y[2].a:=0
    else if k=5
           then delta_y[2].a:=-1
           else begin
                  case k of
                    1 : s:='-0.1';
                    2 : s:='0.8';
                    4 : s:='-0.9';
                  end;
                  delta_y[2].a:=left_read(s);
                end;
  if (k=1) or (k=3)
    then delta_y[2].b:=1
    else begin
           case k of
             2 : s:='1.1';
             4 : s:='0.2';
             5 : s:='-0.7'
           end;
           delta_y[2].b:=right_read(s)
         end;
  if k<=3
    then delta_y[3].a:=-1
    else begin
           case k of
             4 : s:='-0.2';
             5 : s:='0.7'
           end;
           delta_y[3].a:=left_read(s)
         end;
  if k=3
    then delta_y[3].b:=0
    else if k=5
           then delta_y[3].b:=1
           else begin
                  case k of
                    1 : s:='0.1';
                    2 : s:='-0.7';
                    4 : s:='0.9'
                  end;
                  delta_y[3].b:=right_read(s);
                end;
  delta_y[4]:=delta_y[1];
  if k=1
    then Y0[1]:=1
    else begin
           case k of
             2 : s:='0.54030230586806193';
             3 : s:='-0.41614683654775671';
             4 : s:='-0.98999249660561484';
             5 : s:='-0.65364362090705469'
           end;
           Y0[1].a:=left_read(s);
           case k of
             2 : s:='0.54030230586832392';
             3 : s:='0.54030230586832392';
             4 : s:='-0.41614683654635373';
             5 : s:='-0.98999249659487845'
           end;
           Y0[1].b:=right_read(s)
         end;
  if k=1
    then Y0[2]:=0
    else begin
           case k of
             2 : s:='0.84147098480775095';
             3 : s:='0.90929742682451305';
             4 : s:='0.14112000805250249';
             5 : s:='-0.75680249533945081'
           end;
           Y0[2].a:=left_read(s);
           case k of
             2 : s:='0.84147098480799766';
             3 : s:='0.90929742682671995';
             4 : s:='0.14112000806669103';
             5 : s:='-0.75680249527699240'
           end;
           Y0[2].b:=right_read(s)
         end;
  if k=1
    then Y0[3]:=0
    else begin
           case k of
             2 : s:='-0.84147098480813335';
             3 : s:='-0.90929742682673177';
             4 : s:='-0.14112000807149823';
             5 : s:='0.75680249522115091'
           end;
           Y0[3].a:=left_read(s);
           case k of
             2 : s:='-0.84147098480752950';
             3 : s:='-0.90929742682461389';
             4 : s:='-0.14112000804768938';
             5 : s:='0.75680249539553139'
           end;
           Y0[3].b:=right_read(s)
         end;
  if k=1
    then Y0[4]:=1
    else begin
           case k of
             2 : s:='0.54030230586781509';
             3 : s:='-0.41614683654955448';
             4 : s:='-0.98999249661136999';
             5 : s:='-0.65364362092317612'
           end;
           Y0[4].a:=left_read(s);
           case k of
             2 : s:='0.54030230586841365';
             3 : s:='-0.41614683654496195';
             4 : s:='-0.98999249659005160';
             5 : s:='-0.65364362080341119'
           end;
           Y0[4].b:=right_read(s)
         end;
  case k of
    1 : s:='0.99999987500000260';
    2 : s:='0.53988150285540168';
    3 : s:='-0.41660143322387290';
    4 : s:='-0.99006293285765106';
    5 : s:='-0.65326513796975984'
  end;
  Y1[1].a:=left_read(s);
  case k of
    1 : s:='0.99999987500000261';
    2 : s:='0.53988150285566393';
    3 : s:='-0.41660143322246896';
    4 : s:='-0.99006293284690298';
    5 : s:='-0.65326513788159893'
  end;
  Y1[1].b:=right_read(s);
  case k of
    1 : s:='0.00049999997916666';
    2 : s:='0.84174103076555766';
    3 : s:='0.90908923975273205';
    4 : s:='0.14062499418482109';
    5 : s:='-0.75712922253598473'
  end;
  Y1[2].a:=left_read(s);
  case k of
    1 : s:='0.00049999997916667';
    2 : s:='0.84174103076580463';
    3 : s:='0.90908923975494112';
    4 : s:='0.14062499419902004';
    5 : s:='-0.75712922247346684'
  end;
  Y1[2].b:=right_read(s);
  case k of
    1 : s:='-0.00049999997916667';
    2 : s:='-0.84174103076594050';
    3 : s:='-0.90908923975495312';
    4 : s:='-0.14062499420383427';
    5 : s:='0.75712922241758130'
  end;
  Y1[3].a:=left_read(s);
  case k of
    1 : s:='-0.00049999997916666';
    2 : s:='-0.84174103076533594';
    3 : s:='-0.90908923975283270';
    4 : s:='-0.14062499418000074';
    5 : s:='0.75712922259210957'
  end;
  Y1[3].b:=right_read(s);
  if k=1
    then Y1[4]:=Y1[1]
    else begin
           case k of
             2 : s:='0.53988150285515458';
             3 : s:='-0.41660143322567263';
             4 : s:='-0.99006293286340553';
             5 : s:='-0.65326513798591313'
           end;
           Y1[4].a:=left_read(s);
          case k of
            2 : s:='0.53988150285575403';
            3 : s:='-0.41660143322107518';
            4 : s:='-0.99006293284207678';
            5 : s:='-0.65326513786599710'
           end;
           Y1[4].b:=right_read(s)
         end;
  case k of
    1 : s:='0.99999950000004166';
    2 : s:='0.53946056487236853';
    3 : s:='-0.41705592574963295';
    4 : s:='-0.99013312159395930';
    5 : s:='-0.65288649171618417'
  end;
  Y2[1].a:=left_read(s);
  case k of
    1 : s:='0.99999950000004167';
    2 : s:='0.53946056487263114';
    3 : s:='-0.41705592574822784';
    4 : s:='-0.99013312158319907';
    5 : s:='-0.65288649162793551'
  end;
  Y2[1].b:=right_read(s);
  case k of
    1 : s:='0.00099999983333334';
    2 : s:='0.84201086628811103';
    3 : s:='0.90888082540864577';
    4 : s:='0.14012994516089170';
    5 : s:='-0.75745576045021735'
  end;
  Y2[2].a:=left_read(s);
  case k of
    1 : s:='0.00099999983333335';
    2 : s:='0.84201086628835837';
    3 : s:='0.90888082541085727';
    4 : s:='0.14012994517510158';
    5 : s:='-0.75745576038763901'
  end;
  Y2[2].b:=right_read(s);
  case k of
    1 : s:='-0.00099999983333335';
    2 : s:='-0.84201086628849430';
    3 : s:='-0.90888082541086921';
    4 : s:='-0.14012994517992221';
    5 : s:='0.75745576033171060'
  end;
  Y2[3].a:=left_read(s);
  case k of
    1 : s:='-0.00099999983333334';
    2 : s:='-0.84201086628788946';
    3 : s:='-0.90888082540874666';
    4 : s:='-0.14012994515606490';
    5 : s:='0.75745576050638495'
  end;
  Y2[3].b:=right_read(s);
  if k=1
    then Y2[4]:=Y2[1]
    else begin
           case k of
             2 : s:='0.53946056487212131';
             3 : s:='-0.41705592575143449';
             4 : s:='-0.99013312159971258';
             5 : s:='-0.65288649173236853'
           end;
           Y2[4].a:=left_read(s);
           case k of
             2 : s:='0.53946056487272125';
             3 : s:='-0.41705592574683262';
             4 : s:='-0.99013312157837445';
             5 : s:='-0.65288649161230267'
           end;
           Y2[4].b:=right_read(s)
         end;
  Fdelta:=FTY(delta_t,delta_y);
  case k of
    1 : T0:=0;
    2 : T0:=1;
    3 : T0:=2;
    4 : T0:=3;
    5 : T0:=4
  end;
  T1:=T0+h;
  T2:=T1+h;
  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;
  n:=2;
  max_l:=0;
  repeat
    n:=n+1;
    T3:=T2+h;
    T:=T2+hh;
     for i:=1 to 4 do
      Y[i]:=Y2[i]+hh*Fdelta[i];
    i1:=Psi3Explicit(T,Y);
    for i:=1 to 4 do
      begin
        i1[i]:=9*i1[i]-i1[i];
        i2[i]:=(h*h)*(h*h);
        i1[i]:=i2[i]*i1[i]/24
      end;
    i2:=FTY(T2,Y2);
    i3:=FTY(T1,Y1);
    i4:=FTY(T0,Y0);
    for i:=1 to 4 do
      begin
        i2[i]:=7*i2[i]-2*i3[i];
        i2[i]:=i2[i]+i4[i];
        i2[i]:=h*i2[i]/3;
        Y3[i]:=Y1[i]+(i2[i]+i1[i])
      end;
    T:=T3+hh1;
//    iterations
    l:=0;  // number of iterations
    repeat
      l:=l+1;
      find:=True;
      for i:=1 to 4 do
        Y[i]:=Y3[i]+hh1*Fdelta[i];
      i1:=Psi3Implicit(T,Y);
      for i:=1 to 4 do
        begin
          i1[i]:=11*i1[i]-19*i1[i];
          i2[i]:=(h*h)*(h*h);
          i2[i]:=i2[i]*h;
          i1[i]:=i2[i]*i1[i]/720
        end;
      i2:=FTY(T3,Y3);
      i3:=FTY(T2,Y2);
      for i:=1 to 4 do
        i2[i]:=i2[i]+4*i3[i];
      i3:=FTY(T1,Y1);
      for i:=1 to 4 do
        begin
          i2[i]:=i2[i]+i3[i];
          i2[i]:=h*i2[i]/3;
          Y4[i]:=Y1[i]+(i2[i]+i1[i])
        end;
      for i:=1 to 4 do
        if (Y3[i].a=0) or (Y3[i].b=0)
          then if (Abs(Y4[i].a-Y3[i].a)>=1e-18) or (Abs(Y4[i].b-Y3[i].b)>=1e-18)
                 then find:=False  // 1e-18 - accuracy
                 else
          else if (Abs((Y4[i].a-Y3[i].a)/Y3[i].a)>=1e-18)
                   or (Abs((Y4[i].b-Y3[i].b)/Y3[i].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>=2000
             then if next='n'
                    then if n=2000
                           then Results (T3, Y4, 'i')
                           else
                    else begin
                           Results (T3, Y4, 'i');
                           if n<m
                             then Readln
                         end;
    T0:=T1;
    Y0:=Y1;
    T1:=T2;
    Y1:=Y2;
    T2:=T3;
    Y2:=Y3;
    if l>max_l
      then max_l:=l
  until (n=m) or not find;
  if z='f'
    then CloseFile (plik);
  Writeln;
  Writeln ('The end. Press Enter ...');
  Readln
end.