PROGRAM Approximate_Integration
!--------------------------------------------------------------------
! Program to approximate the integral of a function over the interval
! [A,B] using the trapezoidal method.  Identifiers used are:
!   A, B   : the endpoints of the interval of integration
!   N      : the number of subintervals
!   I      : counter
!   DeltaX : the length of the subintervals
!   X      : a point of subdivision 
!   Y      : the value of the function at X
!   Sum    : the approximating sum
!   F      : the integrand (internal function)
!
! Input:  A, B, and N
! Output: Approximation to integral of F on [A, B]
!--------------------------------------------------------------------

  IMPLICIT NONE
  REAL :: A, B, DeltaX, X, Y, Sum
  INTEGER :: N, I

  PRINT *, "Enter the interval endpoints and the # of subintervals:"
  READ *, A, B, N

  ! Calculate subinterval length
  ! and initialize the approximating Sum and X

  DeltaX = (B - A) / REAL(N)
  X = A
  Sum = 0.0

  ! Now calculate and display the sum

  DO I = 1, N - 1
     X = X + DeltaX
     Y = F(X)
     Sum = Sum + Y
  END DO

  Sum = DeltaX * ((F(A) + F(B)) / 2.0 + Sum)

  PRINT '(1X, "Approximate value using", I4, &
        &" subintervals is", F10.5)', N, Sum

CONTAINS

  !-F(X)---------------
  ! The integrand
  !--------------------
  
  FUNCTION F(X)
  
    REAL :: F
    REAL, INTENT(IN) :: X
  
    F = X**2 + 1.0
  
  END FUNCTION F
  
END PROGRAM Approximate_Integration