module integrate interface subroutine open24(f, a, b, integral, error) real, external :: f real, intent(in) :: a, b real, intent(out) :: integral, error end subroutine recursive subroutine adapt(f, a, b, acc, eps, integral, err) real, external :: f real, intent(in) :: a, b, acc, eps real, intent(out) :: integral, err end subroutine end interface end module subroutine open24(f, a, b, integral, error) implicit none real, external :: f real, intent(in) :: a, b real, intent(out) :: integral, error real :: h, f1, f2, f3, f4, integral2 ! Evaluate the function-points: h = b-a f1 = f(a+h*1.0/6.0) f2 = f(a+h*2.0/6.0) f3 = f(a+h*4.0/6.0) f4 = f(a+h*5.0/6.0) ! Calculate the integrals: integral2 = (h/2.0)*(f2+f3) integral = (h/6.0)*(2.0*f1+f2+f3+2.0*f4) write(*,*) integral ! Estimate of the error: error = abs(integral-integral2) end subroutine open24 recursive subroutine adapt(f, a, b, acc, eps, integral, err) use integrate, only : open24 implicit none real, external :: f real, intent(in) :: a, b, acc, eps real, intent(out) :: integral, err real :: tol, integral1, integral2, err1, err2 call open24(f, a, b, integral, err) tol = acc+eps*abs(integral) if (err >= tol) then call adapt(f, a, (a+b)/2.0, acc/sqrt(2.0), eps, integral1, err1) call adapt(f, (a+b)/2.0, b, acc/sqrt(2.0), eps, integral2, err2) integral = integral1+integral2 err = sqrt(err1**2 + err2**2) endif end subroutine adapt