!---------------!
program romberg
implicit none

integer,parameter :: maxstep=10

integer :: i,n,steps
real(8) :: x0,xn,t,polext0,trap(0:maxstep)

write(*,'(a)',advance='no')'Integration bounds x0,xn: '; read*,x0,xn
write(*,'(a)',advance='no')'Number of steps: '; read*,steps
write(*,'(a)',advance='no')'Initial number of intervals: ';read*,n

do i=0,steps-1
  call trapezoidal(x0,xn,n,t,i)
  trap(i)=t
  if (i /= 0) print*,trap(i),polext0(i,trap)
end do

end program romberg
!-------------------!

!-----------------------------------------!
subroutine trapezoidal(x0,xn,n,trap,step)
implicit none

integer :: i,n,step
real(8) :: h,h2,x0,xn,x0h,trap,func

h=(xn-x0)/dble(n*2**step)
if (step == 0) then
  trap=0.5d0*(func(x0)+func(xn))
  do i=1,n-1
    trap=trap+func(x0+h*dble(i))
  end do
else   
  h2=2.d0*h
  x0h=x0+h
  trap=trap/h2
  do i=0,n*2**(step-1)-1
     trap=trap+func(x0h+h2*dble(i))
  end do 
end if
trap=trap*h

end subroutine trapezoidal
!--------------------------!

!--------------------!
function polext0(n,y)
implicit none

integer :: j,k,n
real(8) :: y(0:n),p,polext0

polext0=0.d0
do k=0,n
   p=1.d0
   do j=0,n
      if (j /= k) p=p/(2.d0**(2*(j-k))-1.d0)
   end do
   polext0=polext0+p*y(k)*(-1)**n
end do

end function polext0
!--------------------!

!----------------!
function func(x)
implicit none

real(8) :: x,func

 func=log(x+sqrt(x**2+1))*x**4
! func=x**6

 end function func
!:-----------------!