[交流] [求助]一个关于梯形积分的程序

C234567
      SUBROUTINE TX(FUNC,A,B,N,S)
      IMPLICIT NONE
      INTEGER N, TNM, I
      DOUBLE PRECISION DEL, SUM, X, A, B, S, FUNC
      EXTERNAL FUNC
      IF( N .EQ. 1) THEN
        S = 0.5D0 * (B - A) * (FUNC(A) + FUNC(B))
      ELSE
        TNM = 2**(N - 1)
        DEL = (B - A) / TNM
        X = A
        SUM = 0.0D0
        DO 5 I = 2, TNM
          X = X + DEL
          SUM = SUM + FUNC(X)
5       CONTINUE
        S=0.5D0*(FUNC(A) + FUNC(B) + 2D0 * SUM) * DEL
      ENDIF
      RETURN
      END

      DOUBLE PRECISION FUNCTION FUNC(X)
      IMPLICIT NONE
      DOUBLE PRECISION X
      FUNC = (X**2) * (X**2 - 2.0D0) * SIN(X)
      RETURN
      END

      PROGRAM MAIN
      IMPLICIT NONE
      INTEGER N, I
      DOUBLE PRECISION S, A, B, FUNC
      EXTERNAL FUNC
      A = 0.0D0
      B = 1.5707963D0
      N = 20
      DO 5 I=1,N
        CALL TX(FUNC,A,B,I,S)
        WRITE(*,'(I4,F20.10)') I, S
5     CONTINUE
      END

   1        0.9057727803
   2       -0.0209449905
   3       -0.3614613180
   4       -0.4495837621
   5       -0.4717563216
   6       -0.4773076746
   7       -0.4786960160
   8       -0.4790431327
   9       -0.4791299138
  10       -0.4791516092
  11       -0.4791570331
  12       -0.4791583890
  13       -0.4791587280
  14       -0.4791588128
  15       -0.4791588339
  16       -0.4791588392
  17       -0.4791588406
  18       -0.4791588409
  19       -0.4791588410
  20       -0.4791588410