def update_halley(yvals, y0):
    """Calculate the variable increment using Halley's method.

    Calculate the amount to increment the variable by in one iteration of
    Halley's method. The goal is to find the value of x for which y(x) = y0.
    `yvals` contains the values [y(x0), y'(x0), y''(x0), ...] of the function
    and its derivatives at the current estimate of x. The Halley update
    increment is
        dx0 = (y0 - y)/y'
        dx = dx0 / (1 + dx0*y''/(2*y')).
    Here, dx0 is the increment used by Newton's method. Halley's method
    iterates by adding x += dx at each step.

    Arguments:
        yvals (iterable of float or array): The value of y(x) and its
            derivatives at the current estimate of x. Must contain at least
            3 entries, for y, y' and y''.
        y0 (float or array): The target value of y.

    Returns:
        dx (float or array): The update increment for x.
    """
    dx0 = (y0 - yvals[0])/yvals[1]
    dx = dx0 / (1 + dx0*yvals[2]/(2*yvals[1]))
    return dx