c
c     Numerical Analysis:
c     The Mathematics of Scientific Computing
c     D.R. Kincaid & E.W. Cheney
c     Brooks/Cole Publ., 1990
c
c     Section 7.4
c
c     Computing the integral of a function using Romberg algorithm  
c
c
c     file: romberg.f
c  
      parameter (MM = 4)
      dimension  r(0:MM,0:MM)     
      data a,b/1.0,3.0/
      f(x) = 1.0/x
c
      print *
      print *,' Romberg extrapolation'
      print *,' Section 7.4, Kincaid-Cheney'
      print *
      print 6,'n','R(n,0)','R(n,1)','R(n,2)','R(n,3)','R(n,4)'
c
      h = b-a  
      r(0,0) = 0.5*h*(f(a) + f(b))    
c
      do 4 n = 1,MM
         h = 0.5*h   
         sum = 0.0   
         do 2 i = 1,2**(n-1)      
            sum = sum + f(a + real(2*i-1)*h)   
 2       continue   
         r(n,0) = 0.5*r(n-1,0) + h*sum   
c
         do 3 m = 1,MM
            r(n,m) = r(n,m-1) + (r(n,m-1) - r(n-1,m-1))/((4.**m) - 1.)      
 3       continue
 4    continue
c
      do 5 i=0,MM
         print 7,i,(r(i,j),j = 0,i)
 5    continue
c
 6    format(a6,a13,4(a14))
 7    format(1x,i5,2x,5(e13.6,2x))
      stop
      end