(*  RUNGE-KUTTA-FEHLBERG ALGORITHM 5.3
*
*   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 stepsize HMAX, minimum stepsize HMIN.
*
*  OUTPUT: T, W, H where W approximates Y(T) and stepsize H was exceeded.
*)
TEMP = Input["This is the Runge-Kutta-Fehlberg 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,"Runge-Kutta-Fehlberg Method\n"];
   Write[OUP,"\n"];
   Write[OUP,"T(i)        W(i)        H           R\n"];
   Write[OUP,"\n"]; 
   (* Step 1 *)
   H = HMAX;
   T = A;
   W = ALPHA;
   Write[OUP,SetPrecision[T,9],"      ",SetPrecision[W,9],"          0          0      \n"];
   OK = 1;
   (* Step 2 *)
   While[T < B && OK == 1,
      K1 = H*F[T,W];
      K2 = H*F[T+H/4,W+K1/4];
      K3 = H*F[T+3*H/8,W+(3*K1+9*K2)/32];
      K4 = H*F[T+12*H/13,W+(1932*K1-7200*K2+7296*K3)/2197];
      K5 = H*F[T+H,W+439*K1/216-8*K2+3680*K3/513-845*K4/4104];
      K6 = H*F[T+H/2,W-8*K1/27+2*K2-3544*K3/2565+1859*K4/4104-11*K5/40];
      (* Step 4 *) 
      R = Abs[K1/360-128*K3/4275-2197*K4/75240.0+K5/50+2*K6/55]/H;
      (* Step 5 *) 
      If[R <= TOL,
         (* Step 6 - Approximation accepted *)
         T = T+H;
	 W = W+25*K1/216+1408*K3/2565+2197*K4/4104-K5/5;
         (* Step 7 *)
	 Write[OUP,SetPrecision[T,9],"   ",SetPrecision[W,9]
           ,"      ",SetPrecision[H,9],"       ",SetPrecision[R,9]];
      ];
      (* Step 8 - To avoid underflow *)
      If[R > 1.0 10^-20,
         DELTA = .84*Exp[0.25*Log[TOL/R]],
	 DELTA = 10.0;
      ];
      (* Step 9 - Calculate new H *)
      If[DELTA <= 0.1,
	 H = .1*H,
	 If[DELTA >= 4,
	    H = 4.0*H,
	    H = DELTA*H;
	 ];
      ];
      (* Step 10 *)
      If[H > HMAX,
 	 H = HMAX;
      ];
      (* Step 11 *)
      If[H < HMIN,
         OK = 0,
	 If[T+H > B,
	    If[Abs[B-T] < TOL,
	       T = B,
	       H = B-T;
            ];
         ];
      ];
   ];
   If[OK == 0,
      Write[OUP,"Minimal H exceeded\n"];
   ];
   (* Step 12 - Process is complete *)
   If[OUP == "OutputStream[",NAME," 3]",
      Print["Output file: ",NAME," created successfully\n"];
      Close[OUP]
   ];
];     	       	 	
Quit[];