Improving Newton's iteration where the derivative is near zero? I'm implementing a root-solver for finding x coordinates of a function f(x), after I have an y-coordinate.  The function is periodic, roughly sinusoidal with constant amplitude but non-linearly varying frequency; for an inverse I don't have a closed-form (it is an infinite series), so I use the Newton iteration to find the x-value at a given y beginning the iteration at $x_0$ which is rather near the true value by something like $x=newton(x=x_0,f(x)-y)$.     
In most cases this works fine, however if the y is in the near of a maximum (or minimum) of f, where the shape is very similar to the maximum of a sinus-curve, the newton-iteration does not converge. The wikipedia gives a bit of information about this, but not a workaround. The last way out would be to resort to binary search, but which I'd like to avoid since the computation of f(x) is (relatively) costly.    
Does someone know an improvement in the spirit of the Newton-iteration (which has often quadratic convergence) for this region of y-values?       
[update] hmmm... perhaps it needs only a tiny twist? It just occurs to me, that it might be possible to go along the way how I find the easily approximated x for the maximum y: here I use the Newton on the derivative of the function and search for the zero:
$x_{max}=newton(x=x_0,f(x)') \qquad $ and this has the usual quadratic convergence. But how to apply this for some y in the near of the maximum?
 A: 
In most cases this works fine, however if the y is in the near of a maximum (or minimum) of f, where the shape is very similar to the maximum of a sinus-curve, the newton-iteration does not converge.

It doesn't converge in what way? Knowing the failure mode is important in fixing them. There are three main possibilities: the iteration is converging slowly/asymptotically, the iteration is converging to the wrong value (ie diverging from the answer you want), or the iteration is approaching a limit cycle (bouncing around among 2, 3, or more values repeatedly).
If your problem is that convergence near a minimum is too slow, and the second derivative doesn't vanish there, too, you could try Halley's method. The trouble with Halley's method is that when your state is far from the root it can cause you to overshoot, reducing your basin of convergence. So, I take a cue from "Numerical Recipes in C," when using it I clamp the effect so that it isn't too small or too large:
\begin{align} 
x_{n+1} = x_n - \left(\frac{f_n}{f_n'}\right) \cdot \mathrm{CLAMP}\left(1 - \frac{f_nf_n''}{2[f'_n]^2}, \frac{1}{2}, 2\right)^{-1}.
\end{align}
The ability to speed up convergence by using higher order methods in this class (Househoder's methods) is limited, however, by the number of derivatives that are zero. A simple example of a function that no Newton like method will find the root of quickly: $e^{-1/x^2}$, (all of its derivatives vanish at $x=0$).
If the problem is that the iteration is diverging from the root you want to find then there's a really easy answer: iterate in the other direction. See, in any iterative process like this you will generically have attractors (points the iteration will converged to, if the starting point is in its basin), repulsors (points the iteration will diverge from), and limit cycles (attractors of the double step, triple step, etc). I just recently had a problem with Newton's method always trying to converge to the solution I didn't want in a two solution problem. Flipping the sign to change the direction of iteration fixed that up nicely, because reversed iteration flips the roles of attractors and repulsors.
If the problem is that you're hitting a limit cycle, then that's more tricky. See, if $x_{n+1} = g(x_n)$ is your update rule, then you can define a fixed point as $g(x) = x$. Where the trouble begins is $g$ implies a whole family of update rules, call them $g^{[n]}$, that are the same as applying $g$ $n$ times. The fixed points of $g$ will all be fixed points of $g^{[n]}$, but $g^{[n]}$ will generically have more. Because iterating in $g$ is the same as iterating in $g^{[n]}$, the fixed points of $g^{[n]}$ can become the attractors that dominate the behavior of the sequence, pulling it into a limit cycle of length $n$. If we assume, as is often the case, that the fixed point you want is a reuplsor within the basin of attraction of a limit cycle, then you might be able to find that point if you: 1, detect the cycle; 2, synthesize a guess that is not to close to the limit cycle points (e.g. average them, or use the $x_0$ that converged to the cycle in the first place); and 3, iterate backwards.
I can't offer any guarantees without further information, though, because those basins of attraction are generically fractals in shape.
A: Have you considered implementing a hybrid method? For example, at each step:
IF a Newton step would result in an iteration that is outside the bounds where you have determined the root must lie, then take a bisection step (slower than Newton, but bisection always converges to a root and is not affected by extrema), or a step using a method other than Newton that is not prone to failing near extrema.
ELSE proceed with a Newton step (since it converges quadratically, as you pointed out).
A: If you have some starting point $x_0$, you could obtain $x_1$ by $$x_n=x_{n-1}-\frac{f'(x_{n-1})}{f''(x_{n-1})}$$
and after that use
$$x_n=x_{n-1}-\alpha\frac{f'(x_{n-1})}{f''(x_{n-1})}$$
where $\alpha$ solves $$x_{n-1}-\alpha\frac{f'(x_{n-1})}{f''(x_{n-1})}=x_{n-2}-\alpha\frac{f'(x_{n-2})}{f''(x_{n-2})}$$
For the following function with many $0$-valued derivatives at the minimum $x=1$,
$f(x)=4 (1-x)^8+4 (1-x)^{12}+1888 (1-x)^{14}+1888 (1-x)^{22}$
when starting at $x_0 = 0$ the 3rd iteration overshoots big time, but sufficiently close to the minimum the quadratic convergence begins (see bottom picture, it actually looks super-quadratic) despite 6 derivatives being equal to $0$.
In the case where all derivatives are 0 at the minimum the modification doesn't make the convergence quadratic. I haven't found something that really works in that case.
In higher dimension, where $x_n$ and the fraction are vectors, you can solve the same alpha equation for each parameter, or just those that cause slow convergence.

A: If the minimum/maximum is negative, then an $x_0$ such that $f(x_0)$ is positive is preferred or if the minimum/maximum is positive then an $x_0$ such that $f(x_0)$ negative is more reasonable but if the root is at the minimum/maximum then I don't think there should be a problem.
Something else with a better rate of convergence is the Secant Method. This has a convergence rate of $\cfrac{1+\sqrt 5}{2}=1.618...$ mostly due to the fact that it starts with two initial values.
A: By popular demand from the OP...
In Newton's method you are replacing your function $y=f(x)$ by a linear approximation around the point $x_0$, $y = f(x_0) + f'(x_0) (x-x_0)$, which intersects the x axis ($y=0$) at $x=x_0-f(x_0)/f'(x_0)$.
You could instead approximate by a parabola as $y=f(x_0) + f'(x_0)(x-x_0) +\frac{1}{2}f''(x_0)(x-x_0)^2$, which intercepts the x-axis at $x = x_0 -\frac{f'(x_0)\mp\sqrt{f'(x_0)-2f(x_0)f''(x_0)}}{f''(x_0)}$. You will of course have the issue of having two, not one, possible next iteration point, but there are multiple ways to get around these: choose the closest one, always move up (or down), choose the one with a smallest $f(x_0)$...
A: I couldn't tell you the cost of this idea, but maybe you could work it out:
You are trying to solve for a root of $g$, where $g(x)=f(x)-y$ and you want your solution in a neighborhood of $x_0$. Let's call that solution $x_s$. The problem is that $g'(x)$ is too small near $x_s$, so in the Newton algorithm you divide by very small things, yielding big changes from $x_i$ to $x_{i+1}$; possibly so big that the algorithm converges to a different solution or not at all.  
So what if you engineered a substitute function $\tilde{g}$, who still satisfied the demand that $\tilde{g}(x_s)=0$, but has $\tilde{g}'(x_s)$ not so small. For example, $\tilde{g}(x)$ could equal $g(x)\cdot\ln|g(x)|$. This has $\lim_{x\to x_s}\tilde{g}(x)=0$ and has $\tilde{g}'(x)=g'(x)\cdot(1+\ln|g(x)|)$. The absolute value of $\tilde{g}'(x_s)$ will be quite larger than that of $g'(x_s)$.
A: The issue faced is a case of linear convergence when the root is a higher order root, or at least initially appears this way. This can be handled using an acceleration of Newton's method, commonly given by Aitken's method. The method is as follows:
$$\dot x_n=x_n+\Delta x_n,~\Delta x_n=-\frac{f(x_n)}{f'(x_n)}$$
$$\ddot x_n=\dot x_n+\Delta\dot x_n,~\Delta\dot x_n=-\frac{f(\dot x_n)}{f'(\dot x_n)}$$
$$x_{n+1}=x_n-\frac{(\Delta x_n)^2}{\Delta\dot x_n-\Delta x_n}$$
This will recover most of the lost speed of Newton's method when it struggles. See also this method, which leads to something very similar to MeMyselfI's answer.
If it is reasonable and you expect the root to be a local extrema, you can also resort to optimization techniques instead.
Another approach is to try finding the roots of $f(x)/f'(x)$ i.e. the places where Newton's method will converge 'slowly'. There is, however, the issue that one may have singularities when $f$ is non-zero but $f'$ is, which may make converging to nearby roots difficult. If this works, however, then the full speed of Newton's method is again retained (but at the cost of needing the second derivative).
