C Double integral by two one dimensional subroutines
C
      DOUBLE PRECISION G,ERMAX,XX,EPS,EPS1,A1,B1,ERR,ERR1,RES1
      EXTERNAL G
      COMMON   ERMAX, XX
      INTEGER KF1,IFLAG1
C
      ERMAX = 0.0D0
      EPS   = 1.0D-4
      EPS1  = 9.0D0/10.0D0 * EPS
      A1    = 0.0D0
      B1    = 1.0D0
      CALL DQ1 (G,A1,B1,EPS1,RES1,ERR1,KF1,IFLAG1)
C
      ERR   = ERR1+(B1-A1)*ERMAX
      WRITE (*,'(2X,D25.12,2X,D15.8,4X,2I5)') RES1, ERR, KF1, IFLAG1
      STOP
      END
C
C
      DOUBLE PRECISION FUNCTION G(X)
         DOUBLE PRECISION F,X,ERMAX,XX,EPS2,A2,B2,RES2,ERR2
         EXTERNAL F
         COMMON   ERMAX, XX
C
         XX    = X
         EPS2  = 0.1D0 * 1.D-4
         A2    = 0.0D0
         B2    = 2.0D0
         CALL DQ2 (F,A2,B2,EPS2,RES2,ERR2,KF2,IFLAG2)
C
C      Test IFLAG2
         IF (IFLAG2.NE.0) 
     *        WRITE (*,*) 'ERROR IN Q2, IFLAG, KF2 = ', IFLAG2, KF2
C
         ERMAX = MAX (ERMAX,ERR2)
         G     = RES2
         RETURN
      END
C
C
      DOUBLE PRECISION FUNCTION F(Y)
         DOUBLE PRECISION ERMAX,X,Y
         COMMON   ERMAX, X
C
         F     = EXP (-X*X*Y*Y)
         RETURN
      END