subroutine rombf( a, b, N,MAXN, func, param, R )
integer*4 N, MAXN
real*8 a, b, param(*), R(MAXN,MAXN)
external func
! Function to compute integrals by Romberg algorithm
! R = rombf(a,b,N,MAXN,func,param)
! Inputs
! a,b Lower and upper bound of the integral
! N Romberg table is computed to N by N
! R Array R is dimensioned as R(MAXN,MAXN)
! func Integrand function; the calling sequence
! is: double (*func)( double x, Matrix param )
! param Set of parameters to be passed to function
! Output
! R Romberg table; Entry R(N,N) is best estimate of
! the value of the integral
integer*4 np, i, j, k, m
real*8 h, sumT
!* Compute the first term R(1,1)
h = b - a ! This is the coarsest panel size
np = 1 ! Current number of panels
R(1,1) = h/2 * ( func(a,param) + func(b,param) )
!* Loop over the desired number of rows, i = 2,...,N
do i=2,N
!* Compute the summation in the recursive trapezoidal rule
h = h/2.0 ! Use panels half the previous size
np = 2*np ! Use twice as many panels
sumT = 0.0
do k=1,(np-1),2
sumT = sumT + func( a + k*h, param)
enddo
!* Compute Romberg table entries R(i,1), R(i,2), ..., R(i,i)
R(i,1) = 0.5 * R(i-1,1) + h * sumT
m = 1
do j=2,i
m = 4*m;
R(i,j) = R(i,j-1) + (R(i,j-1) - R(i-1,j-1))/(m-1)
enddo
enddo
return
end