      PROGRAM AREA
************************************************************************
* 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 used                           *
*     I    : counter                                                   *
*     DELX : 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                                             *
*                                                                      *
* Input:  A, B, and N                                                  *
* Output: Approximation to integral of F on [A, B]                     *
************************************************************************

      REAL F, A, B, DELX, 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

      DELX = (B - A) / REAL(N)
      X = A
      SUM = 0.0

* Now calculate and display the sum

      DO 10 I = 1, N - 1
         X = X + DELX
         Y = F(X)
         SUM = SUM + Y
 10   CONTINUE
      SUM = DELX * ((F(A) + F(B)) / 2.0 + SUM)

      PRINT 20, N, SUM
 20   FORMAT (1X, 'APPROXIMATE VALUE USING', I4,
     +            ' SUBINTERVALS IS', F10.5)

      END


** F(X) *********
* The integrand *
*****************

      FUNCTION F(X)
      REAL F, X

      F = X**2 + 1

      END
