Are there special methods for solving $f'(z)=0$ for analytic $f$?

I am trying to solve a bunch of equations for the zeros of the derivative of an analytic function, and I would like to know if there exist methods that exploit this structure to provide better performance than the standard algorithms.

At the moment I am using Mathematica's FindRoot function, which I understand relies on Newton's method. (I iterate over quasirandom seeds as described here.) My problem, however, has additional, structure, of the form $$f'(z)=0$$ for an analytic $f$, so maybe there is a better way to do this. So: are there methods that exploit this structure to provide better performance?

To be a bit more explicit and provide some context in case it's useful: I am writing a Mathematica package to solve the saddle-point equations for high-order harmonic generation (as explained e.g. here), which are of the form \left\{\begin{aligned} (p(t,u)-A(u))^2+\gamma^2 & =0 \\ (p(t,u)-A(t))^2+\gamma^2 & = \omega, \end{aligned}\right. \tag{1} where $p(t,u)=\frac{1}{i\varepsilon+t-u}\int_u^tA(\tau)\mathrm d\tau$ for $\varepsilon$ a small, positive constant, $\gamma>0$ is fixed, and $\omega>0$ is a parameter. Here $A(t)$ is the vector potential of a laser field and can be a three-dimensional vector, but a simple example is $A(t)=A_0\sin(t)$. These equations can be seen as looking for the zeros of both partial derivatives of $$S(t,u)=\gamma^2(t-u)-\omega t+\int_u^t(p(t,u)-A(\tau))^2\mathrm d\tau.$$ $S$ can be assumed to be known explicitly but in some cases it is relatively awkward (with LeafCounts in the >6000 range).

Some notes on the behaviour of this system:

• I am interested in a specific box in the complex plane for $t$ and $u$. In some regimes some roots can wander off the top and bottom of this box (i.e. to higher $|\mathrm{Im}(t)|$ than prescribed), and in those cases the roots would mostly be ignored anyway. It would be nice to have them as it simplifies the accounting, but they're not crucial.
• Some roots which are mostly not of interest can also wander off the sides of the box, and those would definitely get ignored anyway. However, I have no guarantees on the number of roots in my chosen box, and it takes some fiddling after I've got the entire curves w.r.t. $\omega$ to design a function $g(t_s,u_s)$ that will take a root $(t_s,u_s)$ and tell me whether it's a keeper or not.
• The regime $A_0^2\gg\gamma^2$ is of particular interest. In this case the roots will mostly have real $t$ roots, and as $\omega$ varies they will approach each other, have avoided crossings, and veer off into nonzero $\mathrm{Im}(t)$. When $A_0^2\gg\gamma^2$ these avoided crossings can be very tight and happen very quickly with respect to variations in $\omega$, but I'm OK with having to deal with those $\omega$ regions separately.
• At the moment I have an explicit expression for $S(t,u)$ which I differentiate symbolically with Mathematica, and this gives derivatives with LeafCount roughly 6 to 8 times those of $S$. They're not stuff I'd like to simplify by hand.
• I would appreciate methods which integrate well with Mathematica (either already implemented, existent as third-party tools, or relatively easy to implement), but if none are available I'm interested in any methods in this area.
• Sometimes I do have access to reasonable guesses for the roots (via e.g. a related, simpler $A(t)$, or using the results from the previous $\omega$) but I would rather avoid this or at least have methods which work even without such guesses.
• If a simpler model is desired, setting $\gamma=0$ can yield good estimates for $t$ (but not necessarily for $u$).

While it's especially well-suited for polynomial systems, Smale's $\alpha$-theory gives a way to be certain that the Newton iterates for a nonlinear analytic system will converge, and also gives a criterion for when two distinct initial guesses do or do not converge to the same root.
To use it, you'll need both derivatives and bounds on the growth of the Taylor coefficients of the function you're interested in. This is especially nice for polynomials because the Taylor coefficients after the degree are 0. With this in mind, you might be better off transforming your problem by multiplying out the $(i\epsilon + t - u)^{-1}$ term, since that could put a singularity in the complex plane near the region you're interested in. That's just a hunch on my part though.
• Technical note - the $1/(t-u)$ term is actually non-singular, since it multiplies the integral $\int_u^t A(\tau)\mathrm d\tau$ which also vanishes at $t=u$. The small imaginary $i\varepsilon$ is just there to ensure everything plays nice. (Here $p(t,u)$ is actually obtained by taking an explicit $p$ in $S$ and solving the saddle-point equations over $p$ as well, which is much easier than $t$ and $u$ but does have some curious behaviour at small $t-u$.) It is very rare for solutions of $(1)$ to fall in that region. – Emilio Pisanty May 8 '16 at 21:13
• Thanks for the reference to the $\alpha$-theory - I'll give it a careful look. My functions are pretty far from polynomial, though, so I'm wary of asking for things like Taylor coefficients, particularly of high order. Can you comment on what happens in that case? – Emilio Pisanty May 8 '16 at 21:26
• See eqn. (2) in the paper I linked; you need to know $\sup_k\|dg^{-1}d^kg\|^{1/(k-1)}/k!$, where $g = f'$ is the function whose roots you're finding. It's not so much a question of how near/far from polynomial the function is, more of what you know about the decay of the Taylor coefficients. You're right that it could be prohibitively expensive though, so you'll have to see whether this technique is usable to you or not. – Daniel Shapero May 9 '16 at 0:26