faccur Subroutine

    subroutine faccur(me,h0,h1,facc,x0,twoinf,f0,f1)

    implicit none

    class(diff_func),intent(inout) :: me
    real(wp), intent(inout)  :: h0
    real(wp), intent(inout)  :: h1
    real(wp), intent(out)    :: facc
    real(wp), intent(in)     :: x0
    real(wp), intent(in)     :: twoinf
    real(wp), intent(in)     :: f0
    real(wp), intent(in)     :: f1

    real(wp) :: a0,a1,f00,f2,deltaf,t0,t1,df(5)
    integer :: j

    t0 = 0.0_wp
    t1 = 0.0_wp
    if (h0 /= 0.0_wp) then
      if (x0+h0 /= 0.0_wp) then
        f00 = f1
      else
        h0 = 0.875_wp*h0
        f00 = me%f(x0+h0)
      end if
      if (abs(h1) >= 32.0_wp*twoinf) h1 = h1/8.0_wp
      if (16.0_wp*abs(h1) > abs(h0)) h1 = sign(h1,1.0_wp)*abs(h0)/16.0_wp
      if (me%f(x0+h0-h1) == f00) then
        if (256.0_wp*abs(h1) <= abs(h0)) then
          h1 = 2.0_wp*h1
          do
              if (me%f(x0+h0-h1) /= f00 .or. 256.0_wp*abs(h1) > abs(h0)) exit
              h1 = 2.0_wp*h1
          end do
          h1 = 8.0_wp*h1
        else
          h1 = sign(h1,1.0_wp)*abs(h0)/16.0_wp
        end if
      else
        if (256.0_wp*twoinf <= abs(h0)) then
          do
              if (me%f(x0+h0-h1/2.0_wp) == f00 .or. abs(h1) < 4.0_wp*twoinf) exit
              h1 = h1/2.0_wp
          end do
          h1 = 8.0_wp*h1
          if (16.0_wp*abs(h1) > abs(h0)) h1 = sign(h1,1.0_wp)*abs(h0)/16.0_wp
        else
          h1 = sign(h1,1.0_wp)*abs(h0)/16.0_wp
        end if
      end if
    else
      f00 = f0
    end if

    do j = 1,5
      f2 = me%f(x0+h0-real(2*j-1,wp)*h1)
      df(j) = f2 - f00
      t0 = t0+df(j)
      t1 = t1+real(2*j-1,wp)*df(j)
    end do
    a0 = (33.0_wp*t0-5.0_wp*t1)/73.0_wp
    a1 = (-5.0_wp*t0+1.2_wp*t1)/73.0_wp
    facc = abs(a0)
    do j = 1,5
      deltaf = abs(df(j)-(a0+real(2*j-1,wp)*a1))
      if (facc < deltaf) facc = deltaf
    end do
    facc = 2.0_wp*facc

    end subroutine faccur