(*DOUBLE INTEGRAL ALGORITHM 4.4
*
* To approximate I = double integral ( ( f(x,y) dy dx ) ) with
*   limits of integration from a to b for x and from 
*   c(x) to d(x) for y:
*
*    INPUT:   endpoints a, b; positive integers m, n.
*
*   OUTPUT:   approximation J to I.
*)
TEMP=Input["This is Simpson's Method for double integrals.\n
   Input the function F(X,Y) in terms of x and y.\n
   \n
   For example: Cos[x+y]\n"];
F[x_,y_] := Evaluate[TEMP]; 
TEMP=Input["Input the function C(x) in terms of x\n"];
c[x_] := Evaluate[TEMP]; 
TEMP=Input["Input the function D(x) in terms of x\n"];
d[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 a positive integer N\n
   There will be 2N subintervals for the outer integral\n"];
   m=Input["Input a positive integer M\n
   There will be 2M subintervals for the inner integral\n"];
   If[m<=0 || n<=0,
      Input["Integers must be positive\n
      \n
      Press 1 [enter] to continue\n"],
      OK = 1;
   ];
];
If[OK == 1,
   NN=2*n;
   MM=2*m-1;
   (* Step 1 *)
   H=(B-A)/NN;
   (* Use AN, AE, AO, for J(1), J(2), J(3) respectively *)
   (* End terms *)
   AN=0;
   (* Even terms *)
   AE=0;
   (* Odd terms *)
   AO=0;
   (* Step 2 *)
   For[i=0,
      i<=NN,
      i++,
      (* Step 3 - Composite Simpson's Method for X *)
      X=A+i*H;
      YA=c[X];
      YB=d[X];
      HX=(YB-YA)/(2*m);
      (* Use BN, BE, BO for K(1), K(2), K(3) respectively *)
      (* End terms *)
      BN=F[X,YA]+F[X,YB];
      (* Even terms *) 
      BE=0;
      (* Odd terms *)
      BO=0;
      (* Step 4 *)
      For[J=1,
         J<=MM,
	 J++,
         (* Step 5 *)
	 Y=YA+J*HX;
	 Z=F[X,Y];
         (* Step 6 *)
   	 If[Mod[J,2]==0,
	    BE=BE+Z,
	    BO=BO+Z;
	 ];
      ];
      (* Step 7 - use A1 for L, which is the integral of
         F(X(I), Y) from C(X(1)) to D(X(1))
         by Composite Simpson's Method *) 
      A1=(BN+2*BE+4*BO)*HX/3;
      (* Step 8 *)
      If[i==0 || i==NN,
	 AN=AN+A1,
	 If[Mod[i,2]==0,
	    AE=AE+A1,
	    AO=AO+A1;
	 ];
      ];
   ];
   (* Step 9 - Use AC for J *)
   AC=(AN+2*AE+4*AO)*H/3.0;
   (* Step 10 *)
   Print["\n"];
   Print["The integral of F from ",SetPrecision[A,9],
        " to ",SetPrecision[B,9],
      " is ",SetPrecision[AC,9]];
   Print["obtained with N = ",n," and M = ",m];
];
Quit[];