(* SIMPSON'S COMPOSITE ALGORITHM 4.1
*
*  To approximate I=integral ( ( f(x) dx ) ) from a to b:
*
*  INPUT: endpoints a, b; even positive integer n.
*
* OUTPUT: approximation XI to I
*)
TEMP = Input["This is 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,
   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 an even positive integer N.\n"];
   If[ n > 0 && Mod[n,2]==0,
      OK = 1,
      Input["Number must be even and positive\n
      \n
      Press 1 [enter] to continue\n"];
   ];
];
If[OK == 1,
   (* Step 1 *)
   H = (B-A)/n;
   (* Step 2 *)
   Xi0 = F[A]+F[B];
   (* Summation of f(x(2*i-1)) *)
   Xi1 = 0.0; 
   (* Summation of f(x(2*i)) *)
   Xi2 = 0.0; 
   (* Step 3 *)
   NN = n-1;
   For[i = 1,
      i <= NN,
      i++,
      (* Step 4 *)
      X = A+i*H;
      (* Step 5 *)
      If[Mod[i,2] == 0,
         Xi2 = Xi2+F[X],
         Xi1 = Xi1+F[X];
      ];
   ];
   (* Step 6 *)
   Xi = (Xi0+2.0*Xi2+4.0*Xi1)*H/3.0;
   (* Step 7 *)
   Print["\n"];
   Print["The integral of F from ",SetPrecision[A,9]," to ",
         SetPrecision[B,9]
         " is ",SetPrecision[Xi,9]];
];
Quit[];