      SUBROUTINE RADMUL
     1 (F,N,A,B,MINPTS,MAXPTS,EPS,WK,IWK,RESULT,RELERR,NFNEVL,IFAIL)
 
      LOGICAL LDV
 
      DIMENSION A(*),B(*),WK(*)
      DIMENSION CTR(15),WTH(15),WTHL(15),Z(15)
      DIMENSION W(2:15,5),WP(2:15,3)
 
      PARAMETER (R1 = 1, HF = R1/2)
 
      PARAMETER (XL2 =  0.35856 85828 00318 073D0)
      PARAMETER (XL4 =  0.94868 32980 50513 796D0)
      PARAMETER (XL5 =  0.68824 72016 11685 289D0)
 
      PARAMETER (W2 =  980*R1/6561, W4 = 200*R1/19683)
      PARAMETER (WP2 =  245*R1/486, WP4 = 25*R1/729)
 
      DATA (W(N,1),W(N,3),N=2,15)
     1/-0.193872885230909911D+00,  0.518213686937966768D-01,
     2 -0.555606360818980835D+00,  0.314992633236803330D-01,
     3 -0.876695625666819078D+00,  0.111771579535639891D-01,
     4 -0.115714067977442459D+01, -0.914494741655235473D-02,
     5 -0.139694152314179743D+01, -0.294670527866686986D-01,
     6 -0.159609815576893754D+01, -0.497891581567850424D-01,
     7 -0.175461057765584494D+01, -0.701112635269013768D-01,
     8 -0.187247878880251983D+01, -0.904333688970177241D-01,
     9 -0.194970278920896201D+01, -0.110755474267134071D+00,
     A -0.198628257887517146D+01, -0.131077579637250419D+00,
     B -0.198221815780114818D+01, -0.151399685007366752D+00,
     C -0.193750952598689219D+01, -0.171721790377483099D+00,
     D -0.185215668343240347D+01, -0.192043895747599447D+00,
     E -0.172615963013768225D+01, -0.212366001117715794D+00/
 
      DATA (W(N,5),W(N+1,5),N=2,14,2)
     1/ 0.871183254585174982D-01,  0.435591627292587508D-01,
     2  0.217795813646293754D-01,  0.108897906823146873D-01,
     3  0.544489534115734364D-02,  0.272244767057867193D-02,
     4  0.136122383528933596D-02,  0.680611917644667955D-03,
     5  0.340305958822333977D-03,  0.170152979411166995D-03,
     6  0.850764897055834977D-04,  0.425382448527917472D-04,
     7  0.212691224263958736D-04,  0.106345612131979372D-04/
 
      DATA (WP(N,1),WP(N,3),N=2,15)
     1/-0.133196159122085045D+01,  0.445816186556927292D-01,
     2 -0.229218106995884763D+01, -0.240054869684499309D-01,
     3 -0.311522633744855959D+01, -0.925925925925925875D-01,
     4 -0.380109739368998611D+01, -0.161179698216735251D+00,
     5 -0.434979423868312742D+01, -0.229766803840877915D+00,
     6 -0.476131687242798352D+01, -0.298353909465020564D+00,
     7 -0.503566529492455417D+01, -0.366941015089163228D+00,
     8 -0.517283950617283939D+01, -0.435528120713305891D+00,
     9 -0.517283950617283939D+01, -0.504115226337448555D+00,
     A -0.503566529492455417D+01, -0.572702331961591218D+00,
     B -0.476131687242798352D+01, -0.641289437585733882D+00,
     C -0.434979423868312742D+01, -0.709876543209876532D+00,
     D -0.380109739368998611D+01, -0.778463648834019195D+00,
     E -0.311522633744855959D+01, -0.847050754458161859D+00/
 
      RESULT=0
      ABSERR=0
      IFAIL=3
      IF(N .LT. 2 .OR. N .GT. 15) RETURN
      IF(MINPTS .GT. MAXPTS) RETURN
 
      IFNCLS=0
      LDV=.FALSE.
      TWONDM=2**N
      IRGNST=2*N+3
      IRLCLS=2**N+2*N*(N+1)+1
      ISBRGN=IRGNST
      ISBRGS=IRGNST
      IF(MAXPTS .LT. IRLCLS) RETURN
      DO 10 J = 1,N
      CTR(J)=(B(J)+A(J))*HF
   10 WTH(J)=(B(J)-A(J))*HF
 
   20 RGNVOL=TWONDM
      DO 30 J = 1,N
      RGNVOL=RGNVOL*WTH(J)
   30 Z(J)=CTR(J)
      SUM1=F(N,Z)
 
      DIFMAX=0
      SUM2=0
      SUM3=0
      DO 40 J = 1,N
      Z(J)=CTR(J)-XL2*WTH(J)
      F2=F(N,Z)
      Z(J)=CTR(J)+XL2*WTH(J)
      F2=F2+F(N,Z)
      WTHL(J)=XL4*WTH(J)
      Z(J)=CTR(J)-WTHL(J)
      F3=F(N,Z)
      Z(J)=CTR(J)+WTHL(J)
      F3=F3+F(N,Z)
      SUM2=SUM2+F2
      SUM3=SUM3+F3
      DIF=ABS(7*F2-F3-12*SUM1)
      DIFMAX=MAX(DIF,DIFMAX)
      IF(DIFMAX .EQ. DIF) IDVAXN=J
   40 Z(J)=CTR(J)
 
      SUM4=0
      DO 70 J = 2,N
      J1=J-1
      DO 60 K = J,N
      DO 50 L = 1,2
      WTHL(J1)=-WTHL(J1)
      Z(J1)=CTR(J1)+WTHL(J1)
      DO 50 M = 1,2
      WTHL(K)=-WTHL(K)
      Z(K)=CTR(K)+WTHL(K)
   50 SUM4=SUM4+F(N,Z)
   60 Z(K)=CTR(K)
   70 Z(J1)=CTR(J1)
 
      SUM5=0
      DO 80 J = 1,N
      WTHL(J)=-XL5*WTH(J)
   80 Z(J)=CTR(J)+WTHL(J)
   90 SUM5=SUM5+F(N,Z)
      DO 100 J = 1,N
      WTHL(J)=-WTHL(J)
      Z(J)=CTR(J)+WTHL(J)
      IF(WTHL(J) .GT. 0) GO TO 90
  100 CONTINUE
 
      RGNCMP=RGNVOL*(WP(N,1)*SUM1+WP2*SUM2+WP(N,3)*SUM3+WP4*SUM4)
      RGNVAL=W(N,1)*SUM1+W2*SUM2+W(N,3)*SUM3+W4*SUM4+W(N,5)*SUM5
      RGNVAL=RGNVOL*RGNVAL
      RGNERR=ABS(RGNVAL-RGNCMP)
      RESULT=RESULT+RGNVAL
      ABSERR=ABSERR+RGNERR
      IFNCLS=IFNCLS+IRLCLS
 
      IF(LDV) THEN
  110  ISBTMP=2*ISBRGN
       IF(ISBTMP .GT. ISBRGS) GO TO 160
       IF(ISBTMP .LT. ISBRGS) THEN
        ISBTPP=ISBTMP+IRGNST
        IF(WK(ISBTMP) .LT. WK(ISBTPP)) ISBTMP=ISBTPP
       ENDIF
       IF(RGNERR .GE. WK(ISBTMP)) GO TO 160
       DO 130 K = 0,IRGNST-1
  130  WK(ISBRGN-K)=WK(ISBTMP-K)
       ISBRGN=ISBTMP
       GO TO 110
      ENDIF
  140 ISBTMP=(ISBRGN/(2*IRGNST))*IRGNST
      IF(ISBTMP .GE. IRGNST .AND. RGNERR .GT. WK(ISBTMP)) THEN
       DO 150 K = 0,IRGNST-1
  150  WK(ISBRGN-K)=WK(ISBTMP-K)
       ISBRGN=ISBTMP
       GO TO 140
      ENDIF
 
  160 WK(ISBRGN)=RGNERR
      WK(ISBRGN-1)=RGNVAL
      WK(ISBRGN-2)=IDVAXN
      DO 170 J = 1,N
      ISBTMP=ISBRGN-2*J-2
      WK(ISBTMP+1)=CTR(J)
  170 WK(ISBTMP)=WTH(J)
      IF(LDV) THEN
       LDV=.FALSE.
       CTR(IDVAX0)=CTR(IDVAX0)+2*WTH(IDVAX0)
       ISBRGS=ISBRGS+IRGNST
       ISBRGN=ISBRGS
       GO TO 20
      ENDIF
      RELERR=ABSERR/ABS(RESULT)
      IF(ISBRGS+IRGNST .GT. IWK) IFAIL=2
      IF(IFNCLS+2*IRLCLS .GT. MAXPTS) IFAIL=1
      IF(RELERR .LT. EPS .AND. IFNCLS .GE. MINPTS) IFAIL=0
      IF(IFAIL .EQ. 3) THEN
       LDV=.TRUE.
       ISBRGN=IRGNST
       ABSERR=ABSERR-WK(ISBRGN)
       RESULT=RESULT-WK(ISBRGN-1)
       IDVAX0=WK(ISBRGN-2)
       DO 190 J = 1,N
       ISBTMP=ISBRGN-2*J-2
       CTR(J)=WK(ISBTMP+1)
  190  WTH(J)=WK(ISBTMP)
       WTH(IDVAX0)=HF*WTH(IDVAX0)
       CTR(IDVAX0)=CTR(IDVAX0)-WTH(IDVAX0)
       GO TO 20
      ENDIF
      NFNEVL=IFNCLS
      RETURN
      END
