(*  ADAMS VARIABLE STEP-SIZE PREDICTOR-CORRECTOR ALGORITHM 5.5
*
*   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: I, T(I), W(I), H where at the Ith step W(I) approximates.
*)
TEMP = Input["This is the Adams Variable Step-Size Predictor\n
    Corrector Method\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,"ADAMS VARIABLE STEP-SIZE PREDICTOR CORRECTOR METHOD\n"];
   Write[OUP,"\n"];
   Write[OUP,"t           w           h           sigma\n"];
   Write[OUP,"\n"]; 
   (*STEP 2*)
   H = HMAX;
   T[0] = A;
   W[0] = ALPHA;
   (* OK is used in place of FLAG to exit the loop in step 4 *)
   OK = 1;
   (* Done is used in place of last to indicate when last value is calculated *)
   DONE = 0;
   (* STEP 3 *)
   For[KK = 1,
      KK <= 3,
      KK++,
      K1 = H*F[T[KK-1],W[KK-1]];
      K2 = H*F[T[KK-1]+0.5*H,W[KK-1]+0.5*K1];
      K3 = H*F[T[KK-1]+0.5*H,W[KK-1]+0.5*K2];
      K4 = H*F[T[KK-1]+H,W[KK-1]+K3];
      W[KK] = W[KK-1]+(K1+2.0*(K2+K3)+K4)/6.0;
      T[KK] = T[KK-1]+H;
   ];
   (* NFLAG indicates the computation from RK4 *)
   NFLAG = 1;
   i = 5;
   (* Use TT in place of t *)
   TT = T[3]+H;
   (* STEP 4 *)
   While[DONE == 0,
      (* STEP 5 - Predict W(I) *)
      WP = W[i-2]+H*(55.0*F[T[i-2],W[i-2]]-59.0*F[T[i-3],W[i-3]]+37.0*F[T[i-4],W[i-4]]-9.0*F[T[i-5],W[i-5]])/24.0;
      (* Correct W(I) *)
      WC = W[i-2]+H*(9.0*F[TT,WP]+19.0*F[T[i-2],W[i-2]]-5.0*F[T[i-3],W[i-3]]+F[T[i-4],W[i-4]])/24.0;
      SIG = 19.0*Abs[WC-WP]/(270.0*H);
      (* STEP 6 *)
      If[SIG <= TOL,
	 (* STEP 7 - Results accepted *)
         W[i-1] = WC;
	 T[i-1] = TT;
	 (* STEP 8 *)
         If[NFLAG == 1,
	    K=i-3;
	    KK=i-1;
            (* Previous results are also accepted *) 
	    For[J = K,
	       J <= KK,
	       J++,
	Write[OUP,SetPrecision[T[J-1],9],"    ",SetPrecision[W[J-1],9]
                ,"    ",SetPrecision[H,9],"    ",SetPrecision[SIG,9]];
	    ];
	Write[OUP,SetPrecision[T[i-1],9],"    ",SetPrecision[W[i-1],9],
                 "    ",SetPrecision[H,9],"    ",SetPrecision[SIG,9]],
            (* Previous results were already accepted *)
        Write[OUP,SetPrecision[T[i-1],9],"    ",SetPrecision[W[i-1],9],
               "    ",SetPrecision[H,9],"    ",SetPrecision[SIG,9]];
         ];
	 (* STEP 9 *)
         If[OK == 0,
	    DONE = 1,
	    (* STEP 10 *)
            i = i+1;
	    NFLAG = 0;
	    (* STEP 11 *)
            If[SIG <= 0.1*TOL || T[i-2]+H > B,
               (* Increase H if more accuracy than required has been 
                  obtained, or decrease H to include B as a mesh step *)
	       (* STEP 12 - To avoid underflow *)
               If[SIG <= 1.0 10^-20,
		  Q = 4.0,
		  Q = (0.5*TOL/SIG)^(1/4);
	       ];
               (* STEP 13 *)
               If[Q > 4.0,
	   	  H = 4.0*H,
		  H = Q*H;
               ];
	       (* STEP 14 *)
               If[H > HMAX,
		  H = HMAX;
	       ];
	       (* STEP 15 *)
               If[T[i-2]+4.0*H > B,
	 	  H = 0.25*(B-T[i-2]);
		  If[H < TOL,
		     DONE = 1;
		  ];
		  OK = 0;
	       ];
	       (* STEP 16 *)
               For[KK = i-1,
	          KK <= i+2,
		  KK++,
                  K1 = H*F[T[KK-1],W[KK-1]];
                  K2 = H*F[T[KK-1]+0.5*H,W[KK-1]+0.5*K1];
      		  K3 = H*F[T[KK-1]+0.5*H,W[KK-1]+0.5*K2];
     		  K4 = H*F[T[KK-1]+H,W[KK-1]+K3];
     		  W[KK] = W[KK-1]+(K1+2.0*(K2+K3)+K4)/6.0;
     		  T[KK] = T[KK-1]+H;
               ];
	       NFLAG = 1;
	       i = i+3;
            ];
         ],
         (* False branch for step 6 - result rejected *)
         (* STEP 17 *)
         Q = (0.5*TOL/SIG)^(1/4);
	 (* STEP 18 *)
         If[Q < 0.1,
	    H = 0.1*H,
	    H = Q*H;
	 ];
	 (* STEP 19 *)
         If[H < HMIN,
	    Write[OUP,"HMIN exceeded\n"];
	    DONE = 1,
	    If[T[i-2]+4.0*H > B,
	       H = 0.25*(B-T[i-2]);
  	    ];
            If[NFLAG == 1,
	       i = i-3;
 	    ]; 
	    For[KK = i-1,
	       KK <= i+2,
	       KK++, 
               K1 = H*F[T[KK-1],W[KK-1]];
               K2 = H*F[T[KK-1]+0.5*H,W[KK-1]+0.5*K1];
               K3 = H*F[T[KK-1]+0.5*H,W[KK-1]+0.5*K2];
     	       K4 = H*F[T[KK-1]+H,W[KK-1]+K3];
     	       W[KK] = W[KK-1]+(K1+2.0*(K2+K3)+K4)/6.0;
     	       T[KK] = T[KK-1]+H;
            ];
            NFLAG = 1;
            i = i+3;
         ];
      ];
      (* STEP 20 *)
      TT = T[i-2]+H;
   ];
   (* Step 21 *)
   If[OUP == "OutputStream[",NAME," 3]",
      Print["Output file: ",NAME," created successfully\n"];
      Close[OUP]
   ];
];    	 	 	  	 	  	             	  		
Quit[];