code: def find_function_root(fn, x0, xb, delta=0.0, bounds=None):
      
  
      if x0 == xb:
          return np.nan
  
      for i in range(10):
          if np.sign(fn(xb) + delta) != np.sign(fn(x0) + delta):
              break
          if bounds is not None and (xb < bounds[0] or xb > bounds[1]):
              break
          if xb < x0:
              xb *= 0.5
          else:
              xb *= 2.0
  
      # Failed to find a root
      if np.sign(fn(xb) + delta) == np.sign(fn(x0) + delta):
          return np.nan
  
      if x0 == 0:
          xtol = 1e-10 * np.abs(xb)
      else:
          xtol = 1e-10 * np.abs(xb + x0)
  
      return brentq(lambda t: fn(t) + delta, x0, xb, xtol=xtol)
docs: Find the root of a function: f(x)+delta in the interval encompassed
      by x0 and xb.
  
      Parameters
      ----------
  
      fn : function
         Python function.
  
      x0 : float
         Fixed bound for the root search.  This will either be used as
         the lower or upper bound depending on the relative value of xb.
  
      xb : float
         Upper or lower bound for the root search.  If a root is not
         found in the interval [x0,xb]/[xb,x0] this value will be
         increased/decreased until a change in sign is found.