docs: def py_newton(func, x0, fprime=None, args=(), tol=None, maxiter=_iter,
                fprime2=None, low=None, high=None, damping=1.0, ytol=None,
                xtol=1.48e-8, require_eval=False, damping_func=None,
                bisection=False,
                kwargs={}):
      '''Newton's method designed to be mostly compatible with SciPy's 
      implementation, with a few features added and others now implemented.
      
      1) No tracking ofo how many iterations have progressed.
      2) No ability to return a RootResults object
      3) No warnings on some cases of bad input ( low tolerance, no iterations)
      4) Ability to accept True for either fprime or fprime2, which means that
         they are included in the return value of func
      5) No handling for inf or nan!
      
      From scipy, with very minor modifications!
      https://github.com/scipy/scipy/blob/v1.1.0/scipy/optimize/zeros.py#L66-L206
      '''
      if tol is not None:
          xtol = tol
      p0 = 1.0*x0
  #    p1 = p0 = 1.0*x0
  #    fval0 = None
      if bisection:
          a, b = None, None
      if fprime is not None:
          fprime2_included = fprime2 == True
          fprime_included = fprime == True
          
          for iter in range(maxiter):
              if fprime2_included:
                  fval, fder, fder2 = func(p0, *args, **kwargs)
              elif fprime_included:
                  fval, fder = func(p0, *args, **kwargs)
              elif fprime2 is not None:
                  fval = func(p0, *args, **kwargs)
                  fder = fprime(p0, *args, **kwargs)
                  fder2 = fprime2(p0, *args, **kwargs)
              else:
                  fval = func(p0, *args, **kwargs)
                  fder = fprime(p0, *args, **kwargs)
              
              if fder == 0.0 or fval == 0.0:
                  return p0 # Cannot coninue or already finished
              if bisection:
                  if fval < 0.0:
                      a = p0
  #                    f_a = fval
                  else:
                      b = p0
  #                    f_b = fval
              
              fder_inv = 1.0/fder
              # Compute the next point
              step = fval*fder_inv
              if damping_func is not None:
                  if fprime2 is not None:
                      step = step/(1.0 - 0.5*step*fder2*fder_inv)
                  p = damping_func(p0, -step, damping)
  #                variable, derivative, damping_factor
                  
              elif fprime2 is None:
                  p = p0 - step*damping
              else:
                  p = p0 - step/(1.0 - 0.5*step*fder2*fder_inv)*damping
              
              if bisection and a is not None and b is not None:
                  if not (a < p < b) or (b < p < a):
  #                if p < 0.0:
  #                    if p < a:
  #                    print('bisecting')
                      p = 0.5*(a + b)
  #                else:
  #                    if p > b:
  #                        p = 0.5*(a + b)
              
                  
              if low is not None and p < low:
                  p = low
              if high is not None and p > high:
                  p = high
                  
              
  
              
              # p0 is last point (fval at that point), p is new 
                        
              if ytol is not None and xtol is not None:
                  # Meet both tolerance - new value is under ytol, and old value
                  if abs(p - p0) < abs(xtol*p) and abs(fval) < ytol:
                      if require_eval:
                          return p0
                      return p
              elif xtol is not None:
                  if abs(p - p0) < abs(xtol*p):
                      if require_eval:
                          return p0
                      return p
              elif ytol is not None:
                  if abs(fval) < ytol:
                      if require_eval:
                          return p0
                      return p
              
  #            fval0, fval1 = fval, fval0
  #            p0, p1 = p, p0
  # need a history of fval also
              p0 = p
      else:
          return secant(func, x0, args=args, maxiter=maxiter, low=low, high=high,
                        damping=damping,
                        xtol=xtol, ytol=ytol, kwargs=kwargs)
  
      raise UnconvergedError("Failed to converge; maxiter (%d) reached, value=%f " %(maxiter, p))