      PROGRAM ROMBERG
************************************************************************
* Program to approximate the integral of a function over the interval  *
* [A,B] using Romberg integration  Identifiers used are:               *
*     F      : the integrand                                           *
*     A, B   : the endpoints of the interval of integration            *
*     ROM    : two-dimensional array to store Romberg approximations   *
*     LIMIT  : limit on the size of ROM                                *
*     NROWS  : number of rows in the table of Romberg approximations   *
*     NSUBS  : the number of subintervals used                         *
*     K, N   : subscripts                                              *
*     DELX   : the length of the subintervals                          *
*     SUM    : used to obtain one trap. sum from the preceding one     *
*                                                                      *
* Input:  A, B, and NROWS                                              *
* Output: Table of approximations to integral of F on [A, B]           *
************************************************************************

      INTEGER LIMIT, NROWS, NSUBS, N, K
      PARAMETER (LIMIT = 10)
      REAL F, A, B, DELX, SUM, ROM(1:LIMIT, 0:LIMIT)

      PRINT *, 'ENTER THE INTERVAL ENDPOINTS:'
      READ *, A, B
      PRINT *, 'HOW MANY ROWS OF ROMBERG INTEGRATION DO YOU WANT?'
      READ *, NROWS
      PRINT *, '# SUBS                   ROMBERG APPROXIMATIONS'
      PRINT *, '======                   ======================'

* Generate the trapezoidal approximation using 1 subinterval
      DELX = B - A
      ROM(0, 0) = (DELX / 2.0) * (F(A) + F(B))
      NSUBS = 1
      PRINT 10, NSUBS, ROM(0, 0)
 10   FORMAT(1X, I4, 10F12.6)


* Generate Romberg approximations using 2, 4, 8, ... subintervals

      DO 40 N = 1, NROWS - 1

* Generate the trapezoidal approximation using 2**N subintervals
* from that for 2**(N - 1) subintervals

          NSUBS = 2 * NSUBS
          DELX = DELX / 2.0
          SUM = 0.0
          DO 20 K = 1, NSUBS - 1, 2
              SUM = SUM + F(A + K * DELX)
20        CONTINUE
          ROM(N, 0) = ROM(N-1, 0) / 2.0 + DELX * SUM

* Generate the higher-order Romberg approximations for 2**N
* subintervals

          DO 30 K = 1, N
              ROM(N, K) = ROM(N, K-1) +
     +                   (ROM(N, K-1) - ROM(N-1, K-1)) / (4.0**K - 1)
30        CONTINUE

          PRINT 10, NSUBS, (ROM(N, K), K = 1, N)
40    CONTINUE

      END

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

      REAL FUNCTION F(X)
      REAL X

      F = EXP(-X**2)

      END

