(*  ADAMS FOURTH ORDER PREDICTOR-CORRECTOR ALGORITHM 5.4
*
*   To approximate the solution of the initial value problem:
*    Y' = F(T,Y), A <= T <= B, Y(A) = ALPHA,
*    at n+1 equally spaced points in the interval [A,B].
*
*   INPUT: endpoints a, b; initial condition alpha; integer n.
*
*  OUTPUT: approximation W to Y at the (n+1) values of T.
*)
TEMP=Input["This is the Adams-Bashforth Predictor 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,
   n=Input["Input an integer > 3 for the number of\n
   subintervals\n"];
   If[n <= 3,
      Input["Number must be >3. \n;
      \n;
      Press 1 [enter] to continue\n"],
      OK = 1;
   ];
];
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-BASHFORTH FOURTH ORDER PREDICTOR CORRECTOR METHOD\n"];
   Write[OUP,"\n"];
   Write[OUP,"t         w\n"];
   Write[OUP,"\n"]; 
   (* Step 1 *)
   H=(B-A)/n;
   T[0]=A;
   W[0]=ALPHA;
   Write[OUP,SetPrecision[T[0],9],"       ",SetPrecision[W[0],9]];
   (* Step 2 *)
   For[i=1,
      i<=3,
      i++,
      (* Steps 3 and 4 - Compute the starting values using 
         Runge-Kutta method *)
      T[i]=N[T[i-1]+H];
      K1=H*F[T[i-1],W[i-1]];
      K2=H*F[T[i-1]+0.5*H,W[i-1]+0.5*K1]; 
      K3=H*F[T[i-1]+0.5*H,W[i-1]+0.5*K2];
      K4=H*F[T[i],W[i-1]+K3];
      W[i]=W[i-1]+(K1+2.0*(K2+K3)+K4)/6.0;
      (* Step 5 *)
      Write[OUP,SetPrecision[T[i],9],"       ",SetPrecision[W[i],9]];
   ];
   (* Step 6 *)
   For[i=4,
      i <= n,
      i++,
      (* Step 7 - T0 and W0 will be used in place of t and w resp. *)
      T0=N[A+i*H];
      (* Predict W(I) *)
      W0=W[3]+H*(55.0*F[T[3],W[3]]-59.0*F[T[2],W[2]]+37.0*F[T[1],W[1]]-9.0*F[T[0],W[0]])/24.0;
      (* Correct W(I) *)
      W0=W[3]+H*(9.0*F[T0,W0]+19.0*F[T[3],W[3]]-5.0*F[T[2],W[2]]+F[T[1],W[1]])/24.0;
      (* Step 8 *)
      Write[OUP,SetPrecision[T0,9],"       ",SetPrecision[W0,9]];
      (* Step 9 - Prepare for next iteration *)
      For[J=1,
         J<=3,
	 J++,
	 T[J-1]=T[J];
	 W[J-1]=W[J];
      ];
      (* Step 10 *)
      T[3]=T0;
      W[3]=W0;
   ];
];
(* Step 11 *)
If[OUP == "OutputStream[",NAME," 3]",
      Print["Output file: ",NAME," created successfully\n"];
      Close[OUP]
];
Quit[];