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))
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