LR4

      PROGRAM MAIN
      DIMENSION H(10000000)
      COMMON /CONST_PI/ PI
      PI = 3.1415927
      N = 21

C      CALL GRID_CREATION(N, H)
C      CALL GAUSS_3(N, H)

      DO I=1, N

            CALL GRID_CREATION(2**(I-1), H)
C            CALL METHOD_OF_TRAPS(2**(I-1), H)
            CALL SIMPSONS_METHOD(2**(I-1), H)
c            CALL GAUSS_3(2**(I-1), H)
      ENDDO

      END


C ====================================================================
C                             FUNCTION
C ====================================================================

      FUNCTION FUN(X)
      COMMON /CONST_PI/ PI
C      FUN = PI*SIN(2.5*PI*X) + 2.5*PI**2*X*COS(2.5*PI*X)
C      FUN = X**5
      FUN = -(X**5)/120 + (X**3)/6 - X
      RETURN
      END

C ====================================================================
C                           GRID CREATION
C ====================================================================

      SUBROUTINE GRID_CREATION(N, H)
      DIMENSION H(*)

      OPEN(1, FILE = 'input.txt', STATUS = 'OLD')
      READ(1,*) A, B
C      READ(1,*) A, B, N
      CLOSE(1)

      STEP = ((B - A) / N)
      N = N + 1

      DO I = 1, N - 1
            H(I) = A + STEP * (I - 1)
      ENDDO

      H(N) = B
      RETURN
      END


C ====================================================================
C                          METHOD OF TRAPS
C ====================================================================

      SUBROUTINE METHOD_OF_TRAPS(N, H)
      DIMENSION H(*)
      STEP = H(N) - H(N-1)
      PANSW = (FUN(H(1)) + FUN(H(N))) / 2
      DO I = 2, N-1
            PANSW = PANSW + FUN(H(I))
      ENDDO

      PANSW = STEP * PANSW

      PRINT 1, N, PANSW
    1 FORMAT(I5,' Method of traps:',E15.8)
      END

C ====================================================================
C                         SIMPSON'S METHOD
C ====================================================================

      SUBROUTINE SIMPSONS_METHOD(N, H)
      DIMENSION H(*)
      PANSW = 0
      STEP = H(N)-H(N-1)

      DO I=1, N-1
            PANSW = PANSW + STEP * (FUN(H(I)) + 4 * FUN(H(I) + STEP/2)
     +        + FUN(H(I)+STEP)) / 6
      ENDDO

      PRINT 2, N, PANSW
    2 FORMAT(I5,' Simpsons method:', E15.8)
      END

C ====================================================================
C                         GAUSS METHOD
C ====================================================================

      SUBROUTINE GAUSS_3(N, H)
      DIMENSION H(*), W(3), X(3)

      PANSW = 0
      W(1) = 5./9.
      W(2) = 8./9.
      W(3) = 5./9.
      X(1) = -1 * (0.6) ** (0.5)
      X(2) = 0
      X(3) = (0.6) ** (0.5)

      A = H(1)
      B = H(N)
C     DO I = 1, 2
C       PANSW = PANSW + A(I) * DFUN  (1./2 * (B + Aa + X(I)*(B - Aa)))
C     ENDDO
      DO J = 2, N + 1
        DO I = 1, 3
          PANSW = PANSW + W(I) *
     +    FUN(H(J) + 0.5*(X(I)+1)*(H(J-1)-H(J)))
        ENDDO
      ENDDO

      PANSW = PANSW * (B - A) / (2. * N)
      PRINT 3, N, PANSW
    3 FORMAT(I6,' Gauss method:', E15.8)
      END