      subroutine attract_search(xp,yp,zp,root,value,eps)
      double precision xp,yp,zp,root(3),rooti(3),value
      double precision coeff(4),rho,drho(3),d2rho(6),fret
      double precision eps,xpt(3),elfd,delf(3),d2elf(6)
      integer ier,level
      logical func_type
      common/ftype/func_type
      external func
      xpt(1)=xp
      xpt(2)=yp
      xpt(3)=zp
      ier=0
      fret=0.0d0
      if(eps.lt.1.0d0) then
         call dfmfp(func,xpt,fret,eps,ier)
      endif
      xp=xpt(1)
      yp=xpt(2)
      zp=xpt(3)
      level=2
      call elf_eval(rho,drho,d2rho,elfd,delf,d2elf,xp,yp,zp,level)
      if(func_type) then
         coeff(1)=1.0d0
         coeff(2)=-(d2elf(1)+d2elf(4)+d2elf(6))
         coeff(3)=d2elf(1)*d2elf(4)+d2elf(1)*d2elf(6)+
     .      d2elf(4)*d2elf(6)-d2elf(2)*d2elf(2)-
     .      d2elf(3)*d2elf(3)-d2elf(5)*d2elf(5)
         coeff(4)=d2elf(1)*d2elf(5)*d2elf(5)+
     .      d2elf(4)*d2elf(3)*d2elf(3)+
     .      d2elf(6)*d2elf(2)*d2elf(2)-
     .      d2elf(1)*d2elf(4)*d2elf(6)-
     .      2.0d0*d2elf(2)*d2elf(3)*d2elf(5)
         value=elfd
      else
         coeff(1)=1.0d0
         coeff(2)=-(d2rho(1)+d2rho(4)+d2rho(6))
         coeff(3)=d2rho(1)*d2rho(4)+d2rho(1)*d2rho(6)+
     .      d2rho(4)*d2rho(6)-d2rho(2)*d2rho(2)-
     .      d2rho(3)*d2rho(3)-d2rho(5)*d2rho(5)
         coeff(4)=d2rho(1)*d2rho(5)*d2rho(5)+
     .      d2rho(4)*d2rho(3)*d2rho(3)+
     .      d2rho(6)*d2rho(2)*d2rho(2)-
     .      d2rho(1)*d2rho(4)*d2rho(6)-
     .      2.0d0*d2rho(2)*d2rho(3)*d2rho(5)
         value=rho
      endif
      call cubic(coeff,root,rooti)
      return
      end
