(* ROMBERG ALGORITHM 4.2
*
*  INPUT: endpoints a, b; integer n. 
*
* OUTPUT: an array R.  ( R(2,n) is the approximation to I. )
*     R is computed by rows;only two rows saved in storage
*)
TEMP=Input["This is Romberg Integration.\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,
   A = Input["Input the lower limit of integration\n"];
   B = Input["Input the upper limit of integration\n"];
   If[A > B,
      Input["Lower limit must be less than upper limit\n
      \n
      Press 1 [enter] to continue\n"],
      OK = 1;
   ];
];
OK = 0;
While[OK == 0,
   n = Input["Input the number rows - no decimal point\n"];
   If[n > 0,
      OK = 1,
      Input["Number must be positive\n
      \n
      Press 1 [enter] to continue\n"];
   ];
];
(* Step 1*)
If[OK == 1,
   H = B-A;
   R[0,0] = (F[A]+F[B])/2*H;
   (* Step 2 *)
   Print["Initial data:\n"];
   Print["Limits of integration = [ ",A," , ",B," ]"];
   Print["Number of rows = ",n];
   Print["\n"];
   Print["Romberg Integration Table:\n"];
   Print["\n"];
   Print[SetPrecision[R[0,0],9]];
   Print["\n"];
   (* Step 3 *)
   For[i = 2,
      i <= n,
      i++,
      (* Step 4 - Approximation from Trapezoidal method *)
      TOTAL = 0;
      M = 2^(i-2);
      For[K = 0,
         K <= M-1,
         K++;
        TOTAL = TOTAL+F[A+(K-0.5)*H];
      ];
      R[1,0] = (R[0,0]+H*TOTAL)/2;
      (* Step 5 - Extrapolation *)
      For[J = 2,
         J <= i,
	 J++,
   	 L = 2^(2*(J-1));
     	 R[1,J-1] = R[1,J-2]+(R[1,J-2]-R[0,J-2])/(L-1);
      ];
      (* Step 6 *)
      For[K = 1,
         K <= i,
	 K++,
	 Print[SetPrecision[R[1,K-1],9]];
      ];
      Print["\n"];
      Print["\n"];
      (* Step 7 *)
      H = H/2;
      (* Since only two rows are kept in storage, this step
         is to prepare for the next row 1 to I of R 
         Step 8 *)
      For[J = 1,
         J <= i,
	 J++,
	 R[0,J-1] = R[1,J-1];
      ];
   ];
];
Quit[];