(*  EXTRAPOLATION ALGORITHM 5.6
*
*   To approximate the solution of the initial value problem:
*    Y' = F(T,Y), A <= T <= B, Y(A) = ALPHA,
*    with local truncation error within a given tolerance.
*
*   INPUT: endpoints a, b; initial condition alpha; tolerance TOL;
*	    maximum step size HMAX; minimum step size HMIN.
*
*  OUTPUT: T, W, H where W approximates y(T) and stepsize H was
*           used or a message that the minimum stepsize was exceeded.
*)
TEMP=Input["This is Gragg Extrapolation\n
  Input the function F(T,Y) in terms of t and y\n 
 \n
 For example: y-t^2+1\n"];
F[t_,y_] :=Evaluate[TEMP];
OK = 0;
While[OK == 0,
   A=Input["Input the left endpoint\n"];
   B=Input["Input the right endpoint\n"];
   If[A >= B,
      Input["Left endpoint must be less than right endpoint\n
      \n
      Press 1 [enter] to continue\n"],
      OK = 1;
   ];
];
ALPHA=Input["Input the initial condition\n"];
OK=0;
While[OK == 0,
   TOL=Input["Input tolerance\n"];
   If[TOL <= 0,
      Input["Tolerance must be positive\n
      \n
      Press 1 [enter] to continue\n"],
      OK = 1;
   ];
];
OK = 0;
While[OK == 0,
   HMIN=Input["Input minimum mesh spacing\n"];
   HMAX=Input["Input maximum mesh spacing\n"];
   If[HMIN < HMAX && HMIN > 0,
      OK = 1,
      Input["Minimum mesh spacing must be a positive real number\n
      and less than the maximum mesh spacing\n
      \n
      Press 1 [enter] to continue\n"];
   ];
];
If[OK == 1,
   FLAG = Input["Select output destination\n
                 1. Screen\n
                 2. Text file\n
                 Enter 1 or 2\n"];
   If[FLAG == 2,
      NAME = InputString["Input the file name\n
                          For example:   output.dta\n"];
      OUP = OpenWrite[NAME,FormatType->OutputForm],
      OUP = "stdout";
   ];
   Write[OUP,"GRAGG EXTRAPOLATION\n"];
   Write[OUP,"\n"];
   Write[OUP,"T           W          H           K\n"];
   Write[OUP,"\n"]; 
   (* Step 1 *)
   NK[0]=2;
   NK[1]=4;
   For[J=1,
      J <= 3,
      J++,
      i=2*J;
      NK[i]=3*NK[i-1]/2;
      NK[i+1]=2*NK[i-1];
   ];
   (* Step 2 *)
   H=HMAX;
   T0=A;
   W0=ALPHA;
   (* Done is used in place of FLAG to exit loop in Step 4 *)
   DONE = 0;
   (* Step 3 *)
   For[i=1,
      i <= 7,
      i++,
      For[J = 1,
         J <= i,
	 J++,
	 q[i-1,J-1]=(NK[i]*1/NK[J-1])*(NK[i]*1/NK[J-1]);
      ];
   ];
   (* Step 4 *)
   While[DONE == 0,
      (* Step 5 *)
      K=1;
      (* When desired accuracy achieved, NFLAG is set to 1 *)
      NFLAG = 0;
      (* Step 6 *)
      While[K <= 8 && NFLAG == 0,
         (* Step 7 *)
         HK=H/NK[K-1];
   	 T=T0;
	 W2=W0;
         (* Euler first step *)
	 W3=W2+HK*F[T,W2];
	 T=T0+HK;
         (* Step 8 *)
	 M=NK[K-1]-1;
         For[J=1,
	    J <= M,
	    J++,
	    W1=W2;
	    W2=W3;
            (* Midpoint method *)
	    W3=W1+2*HK*F[T,W2];
	    T=T0+(J+1)*HK;
         ];
         (* Step 9 - Endpoint correction to compute Y(K, 1) *)
         Y[K]=(W3+W2+HK*F[T,W3])/2;
         (* Step 10 *)
         (* Note : Y(K-1)=Y(K-1,1), Y(K-2)=Y(K-1,2),...,
            Y(1)=Y(K-1,K-1) since only previous row of table is saved *)
	 If[K >= 2,
            (* Step 11 *)
	    J=K;
            (* Save Y(K-1,K-1) *) 
	    V=Y[1];
            (* Step 12 *)
	    While[J >= 2,
               (* Extrapolation to compute Y(J-1)=Y(K,K-J+2) *) 
	       Y[J-1]=Y[J]+(Y[J]-Y[J-1])/(q[K-2,J-2]-1);
	       J=J-1;
	    ];
            (* Step 13 *)
	    If[Abs[Y[1]-V] <= TOL,
	       NFLAG = 1;
	    ];
         (* Y(1) accepted as new w *)
	 ];
         (* Step 14 *)
	 K=K+1;
      ];
      (* Step 15 *)
      K=K-1;
      (* Step 16 *)
      If[NFLAG == 0,
         (* Step 17 - New value for w rejected, decrease H *)
	 H=H/2;
         (* Step 18 *)
	 If[H < HMIN,
	    Write[OUP,"HMIN exceeded\n"];
	    DONE = 1;
	 ],
         (* Step 19 - New value for w accepted *) 
	 W0=Y[1];
	 T0=T0+H;
	 Write[OUP,SetPrecision[T0,9],"         ",SetPrecision[W0,9]
         ,"        ",SetPrecision[H,9],"          ",K];
         (* Step 20 - Increase H if possible *)
         If[T0 >= B,
            DONE = 1,
            If[T0+H > B,
               H= B-T0,
	       If[K <= 3,
                  H=2*H;
	       ];
	    ];
         ];
	 If[H > HMAX,
            H=H/2;
         ];
      ];
   ];
   (* Step 21 *)
   If[OUP == "OutputStream[",NAME," 3]",
      Print["Output file: ",NAME," created successfully\n"];
      Close[OUP];
   ];
];
Quit[];