(* ADAPTIVE QUADRATURE ALGORITHM 4.3
*
*  To approximate I=integral ( ( f(x) dx ) ) from a 
*    to b within a given tolerance TOL:
*
*  INPUT: endpoints a, b; tolerance TOL; 
*    limit N to number of levels
*
* OUTPUT: approximation APP or message that N is exceeded
*)
TEMP = Input["This is Adaptive Quadrature with Simpson's Method.\n
   Input the function F(X) in terms of x.\n
   \n
   For example: Cos[x]\n"];
F[x_] := Evaluate[TEMP];
OK = 0;
While[OK == 0,
   AA = Input["Input the lower limit of integration\n"];
   BB = Input["Input the upper limit of integration\n"];
   If[AA > BB,
      Input["Lower limit must be less than upper limit\n
      \n
      Press 1 [enter] to continue\n"],
      OK = 1;
   ];
];
OK = 0;
While[OK == 0,
   EPS = Input["Input tolerance\n"];
   If[EPS > 0,
      OK = 1,
      Input["Tolerance must be positive\n
      \n
      Press 1 [enter] to continue\n"];
   ];
];
OK = 0;
While[OK == 0,
   n = Input["Input the maximum number of levels\n"];
   If[n > 0,
      OK = 1,
      Input["Number must be positive\n
      \n
      Press 1 [enter] to continue\n"];
   ];
];
If[OK == 1,
   CNT = 0;
   OK = 1;
   (* Step 1 *)
   APP = 0;
   i = 1;
   TOL[i] = 10*EPS;
   A[i] = AA;
   H[i] = 0.5*(BB-AA);
   FA[i] = N[F[AA]];
   CNT = CNT+1;
   FC[i] = N[F[(AA + H[i])]];
   CNT = CNT+1;
   FB[i] = N[F[BB]];
   CNT = CNT+1;
   (* Approximation from Simpson's method for entire interval *)
   S[i] = H[i]*(FA[i]+4*FC[i]+FB[i])/3.0;
   L[i] = 1;
   (* Step 2 *)
   While[i > 0 && OK == 1,
      (* Step 3 *)
      FD = N[F[(A[i]+0.5*H[i])]];
      CNT = CNT+1; 
      FE = N[F[(A[i]+1.5*H[i])]];
      CNT = CNT+1;
      (* Approximations from Simpson's method for halves of intervals *)
      S1 = H[i]*(FA[i]+4*FD+FC[i])/6.0;
      S2 = H[i]*(FC[i]+4*FE+FB[i])/6.0;
      (* Save data at this level *)
      V[0] = A[i];
      V[1] = FA[i];
      V[2] = FC[i];
      V[3] = FB[i];
      V[4] = H[i];
      V[5] = TOL[i];
      V[6] = S[i];
      LEV = L[i];
      (* Step 4 - Delete the level *)
      i = i-1;
      (* Step 5 *)
      If[Abs[S1+S2-V[6]] < V[5],
	 APP = APP+(S1+S2),
	 If[LEV >= n,
	    OK = 0,
	    (* Procedure fails - Add one level *)
	    (* Data for right half of subinterval *)
	    i = i+1;
	    A[i] = V[0]+V[4];
	    FA[i] = V[2];
	    FC[i] = FE;
	    FB[i] = V[3];
	    H[i] = 0.5*V[4];
	    TOL[i] = 0.5*V[5];
	    S[i] = S2;
	    L[i] = LEV+1;
	    (* Data for left half subinterval *)
	    i = i+1;
	    A[i] = V[0];
	    FA[i] = V[1];
	    FC[i] = FD;
	    FB[i] = V[2];
	    H[i] = H[i-1];
	    TOL[i] = TOL[i-1];
	    S[i] = S1;
	    L[i] = L[i-1];                     
	 ];
      ];
   ];
   If[OK == 0,
      Print["Level exceeded.  Method did not produce an \n"];
      Print["accurate approximation.\n"],
      Print["\n"];
      Print["The integral of F from ",SetPrecision[AA,9],
           " to ",SetPrecision[BB,9]," is"];
      Print[SetPrecision[APP,9],"  to within \n",EPS];
      Print["The number of function evaluations is: \n",CNT];
   ];
];
Quit[];