各位好，我用变步长sipson子例行程序计算一个二重积分，结果总是不对。如果换别的被积函数，算出的结果没问题。说明这个二重积分的子程序时没错的。但是我就是找不到问题出在哪里。请大家帮忙看看好吗？下面是代码：
    PROGRAM MAIN
    COMMON /AB1/A
    COMMON/AB2/B
    COMMON/R/R
    PARAMETER(MAXR=120)
    REAL ::RESULT1
    EXTERNAL RT
    PARAMETER(THETA=0.523599)
    PARAMETER(MASS=131)
    PARAMETER(PI=3.1415926)
    
    B=RT(0.0)
    A=RT(1.5707963)
    EPS=1.E-4

    GS=1./(SQRT(2*PI))**3
    FA=(3*MASS)/(2*B*(A**2))
    FAGS=FA*GS

    write(*,*)'FAGS=',FAGS
    WRITE(*,*)'A=',A,'B=',B
    write(*,*)'mass=',mass


       CALL SIMPS2(0.0,PI,EPS,KFL,NNN,RESULT1)
        WRITE(*,*) R,RESULT1*FAGS
     

    END
***********************
    BLOCK DATA
    COMMON /R/R
        DATA R/1./
    END BLOCK DATA
************************
    SUBROUTINE FLS(THETA1,Y1,Y2)
    COMMON/AB1/A
    COMMON/AB2/B
    EE=(A**2-B**2)/(B**2)
    Y1=0.
    Y2=A*(1-0.5*EE*(COS(THETA1)**2)+0.75*(EE**2)*(COS(THETA1)**4))
    END SUBROUTINE
************************************
    REAL FUNCTION F(R1,THETA1)
    PARAMETER(THETA=1.5707963)
    PARAMETER(MASS=131)
    COMMON/R/R
    EXTERNAL BESSI0
    
    F=EXP(-0.5*(R1**2+R**2))*EXP(R*R1*COS(THETA)*COS(THETA1))*  
     &    SIN(THETA1)*(R1**2)*

    END FUNCTION
*************************************
    REAL FUNCTION  RT(THE)
    PARAMETER(THETA=1.5707963)
    PARAMETER(MASS=131)
    PARAMETER(PI=3.1415926)
    BETA2=-0.113    
    AA=0.523
    Y20=0.25*SQRT(5./PI)*(-1+3*COS(THE)**2.)        
    SR0=(1.23*MASS**(1./3.)-0.6)**2
    R0=SQRT(SR0+(7./3.)*(AA**2)*(PI**2)-5*(0.9**2))
    RT=R0*(1+BETA2*Y20)
    END FUNCTION
******************************************************

***************************************
        SUBROUTINE SIMPS1(X,EPS,KFL,N1,S)

    N1=1
    CALL FLS(X,Y1,Y2) 
    H1=0.5*(Y2-Y1)
    E=ABS(Y1)+ABS(Y2)
    T1=H1*(F(X,Y1)+F(X,Y2))
1    Y=Y1-H1
    T2=0.5*T1
     DO 2 I=1,N1
    Y=Y+2.*H1
2    T2=T2+H1*F(X,Y)
    S=(4.0*T2-T1)/3.0
    N1=N1+N1
    IF(N1.LT.8)GOTO 3
    IF(ABS(S-S0).LT.(1.0+ABS(S))*EPS)RETURN
3    S0=S
    T1=T2
    H1=0.5*H1
    IF(E+H1.NE.E)GOTO 1
    IF(Y1.EQ.Y2)RETURN
    KFL=KFL+1
    RETURN
    END

    
************************************************************
    SUBROUTINE SIMPS2(A,B,EPS,KFL,N,SUM)

    N2=1
    KFL=0
    C=ABS(A)+ABS(B)    
    H2=0.5*(B-A)
    CALL SIMPS1(A,EPS,KFL,N1,SS1)

    CALL SIMPS1(B,EPS,KFL,N3,SS2)

    N=N1+N3+2
    TB1=H2*(SS1+SS2)
4    X=A-H2
    TB2=0.5*TB1
    DO 5 I=1,N2
       X=X+2.0*H2
    CALL SIMPS1(X,EPS,KFL,N1,S)
    N=N+N1+1
5    TB2=TB2+H2*S
    SUM=(4.0*TB2-TB1)/3.0
    N2=N2+N2
    IF(N2.LT.8)GOTO 6
    IF(ABS(SUM-SUM0).LE.(1.0+ABS(SUM))*EPS) RETURN
6    SUM0=SUM
    TB1=TB2
    H2=0.5*H2
    IF(C+H2.NE.C)GOTO 4
    N=-N
    RETURN
    END
这个程序算出的结果是0.0659，但是我用mathematica算出的结果是2.21208，后者的结果应该是正确的。实在找不出原因。请大家帮忙看看好吗？多谢了