bisect Subroutine

subroutine bisect(f,x1,x2,eps,niter,flag)
  implicit none
  interface                                                             
     function f(x)                                                   
       real(8) :: x
       real(8) :: f
     end function f
  end interface
  real(8)          :: x1, x2
  real(8),optional :: eps
  integer,optional :: Niter
  integer,optional :: flag
  real(8)          :: Root
  real(8)          :: eps_
  integer          :: Niter_
  real(8)          :: a, b, c
  integer          :: iter
  !
  eps_=1d-9;if(present(eps))eps_=eps
  Niter_=200;if(present(Niter))Niter_=Niter
  if(f(x1)*f(x2)>0.0) then
     flag = 0
     return
  end if
  a=x1
  b=x2
  do iter=1,Niter_
     c=(b+a)/2d0
     if(f(a)*f(c) <= 0d0) then
        b = c
     else
        a = c
     end if
     if(abs(b-a) <= eps_) exit  
  end do
  root=(b+a)/2.0
  !check if it is a root or singularity
  if(present(flag))then
     if (abs(f(root)) < 1d0) then
        flag =  iter
     else
        flag = -iter
     end if
  endif
  x1=root
  x2=root
end subroutine bisect