Here is my approach using an explicit pseudo-arclength method.
The typical predictor-corrector method uses the parameter as continuation parameter. For example, let
$$
G(u,h) = 0. \label{eq:root}\tag{1}
$$
If we know a solution $(u_0,h_0)$, then we predict the solution at $h_1 = h_0 + \Delta h$ by noting that, if the Jacobian matrix $G_u$ is invertible, then
$$
u'(h_0) = - G_u^{-1}(u(h_0),h_0)G_h(u(h_0),h_0) = u_0'.
$$
Then we propose the predictor
$$
u_p = u_0 + u_0' \Delta h.
$$
and use some method (usually Newton's) to obtain the corrector.
The problem with this method is that it fails in a fold. To circumvent this problem, instead of parameterizing the solution curve by the continuation parameter, we do it by its arclength. In this case, you have to think that everything is a function of $s$, i.e., $X(s) = \big(u(s),h(s)\big)$, like so (taken from [1]):
In this scheme, not only $u$ is unknown, but also $h$. There are several ways to close the system, the most common being to use the scalar normalization
$$
(X_1 - X_0)^T \dot{X_0} = \Delta s.
$$
This is the equation of a plane, which is perpendicular to the tangent $\dot{X}(s)$ at a distance $\Delta s$ from $X_0$. This plane will intersect the solution curve if $\Delta s$ and the curvature of the curve is not to large. So, we extend \eqref{eq:root} to
$$
G(u,h) = 0,
$$
$$
(u - u_0)^T \dot{u}_0 + (h - h_0)\dot{h}_0 - \Delta s = 0,
$$
where $(u_0,h_0) = \big(u(s_0),h(s_0)\big)$ is a known solution. Now we only need to calculate $\dot{u}_0$ and $\dot{h}_0$ to get things going. This can be done by differentiating $G$ with respect to $s$ and using the normalization condition $\dot{u}_0^T \dot{u}_0 + \dot{h}_0^2 = 1$:
\begin{align}\label{eq:extroot}\tag{2}
\dot{u}_0 &= -G_u^{-1}(u_0,h_0) G_h(u_0,h_0) \dot{h}_0,\\
\dot{h}_0 &= \pm\left(1+\|G_u^{-1}(u_0,h_0) G_h(u_0,h_0)\|^2\right)^{-1/2},
\end{align}
where we choose the sign depending on the direction we want to go. In this scheme, given a solution $(u_0,h_0)$, our predictor is
\begin{align}
u_p &= u_0 + \dot{u}_0 \Delta s,\\
h_p &= h_0 + \dot{h}_0 \Delta s.
\end{align}
Finally, a practical way to approximate the derivatives of $u$ and $h$ is, instead of solving \eqref{eq:extroot} at each step and determining the sign, approximate this equation by
$$
\pmatrix{G_u(u_1,h_1) & G_h(u_1,h_1) \\ \dot{u}_0 & \dot{h}_0}\pmatrix{\dot{u}_1 \\ \dot{h}_1} = \pmatrix{0 \\ 1},
$$
which has the advantage of choosing the right direction at each step.
Details on why this scheme works on folds can be found in Keller's classic notes Lectures on Numerical Methods In Bifurcation Problems.
Some Mathematica examples:
I'm not an MMA expert; these codes are for illustrative purposes only, and haven't been crafted for speed or elegance. I'd love to see some of our house names take a shot into packing TrackRoot[..., Method -> PseudoArcLength] ;)
OP's own example:
funG = {x^2 + y^2 - h, x y - 24};
extfunG = {x^2 + y^2 - h, x y - 24, (x - xPrev) dxPrev + (y - yPrev) dyPrev +
(h - hPrev) dhPrev - ds};
chi = LinearSolve[D[funG, {{x, y}}], -D[funG, h]];
extchi = LinearSolve[Join[D[funG, {{x, y, h}}], {{dxPrev, dyPrev, dhPrev}}], {0, 0, 1}];
{x0, y0, h0} = {x, y, 100} /. FindRoot[funG /. h -> 100, {x, 4}, {y, 9}, MaxIterations->500];
dh0 = -(1 + chi.chi)^(-1/2) /. {x -> x0, y -> y0, h -> h0}; (* left continuation *)
{dx0, dy0} = -chi dh0 /. {x -> x0, y -> y0, h -> h0};
ClearAll[predcorr];
predcorr[xi_, dxi_, step_] := Module[{xp, yp, hp},
{xp, yp, hp} = xi + step dxi;
{{x, y, h}, extchi /. {dxPrev->dxi[[1]], dyPrev->dxi[[2]], dhPrev->dxi[[3]]}, step} /.
FindRoot[extfunG /. {xPrev -> xi[[1]], yPrev -> xi[[2]], hPrev -> xi[[3]],
dxPrev -> dxi[[1]], dyPrev -> dxi[[2]], dhPrev -> dxi[[3]], ds -> step},
{x, xp}, {y, yp}, {h, hp}, MaxIterations -> 500]
]
We track the root using a NestWhileList. Note the Check function, ready to diminish the step-size in order to ensure convergence (There might be a more elegant way to address this problem).
solCurve = NestWhileList[
Check[predcorr[#1, #2, #3], {#1, #2, #3/2}] & @@ # &,
{{x0, y0, h0}, {dx0, dy0, dh0}, 0.1}, #[[1, 3]] < 101 &
] // DeleteDuplicates[#, (#1[[1]] == #1[[2]] &)] &;
Here is a bifurcation diagram:
ListPlot[{#1[[3]], #1[[1]]} & @@@ solCurve, AxesLabel -> {h, x}]
Chris K example:
funG1 = -z^3 + 5/2 z + h;
extfunG1 = {-z^3 + 5/2 z + h, (z - zPrev) dzPrev + (h - hPrev) dhPrev - ds};
chi1 = -D[funG1, h]/D[funG1, z];
extchi1 = LinearSolve[Join[{D[funG1, {{z, h}}]}, {{dzPrev, dhPrev}}], {0, 1}];
{z0, h0} = {z, -10} /. FindRoot[funG1 /. h -> -10, {z, -3}, MaxIterations -> 500];
dh0 = (1 + chi1^2)^(-1/2) /. {z -> z0, h -> h0}; (* right continuation *)
dz0 = -chi1 dh0 /. {z -> z0, h -> h0};
ClearAll[predcorr1];
predcorr1[zi_, dzi_, step_] := Module[{zp, hp},
{zp, hp} = zi + step dzi;
{{z, h}, extchi1 /. {dzPrev -> dzi[[1]], dhPrev -> dzi[[2]]}, step} /.
FindRoot[extfunG1 /. {zPrev -> zi[[1]], hPrev -> zi[[2]], dzPrev -> dzi[[1]],
dhPrev -> dzi[[2]], ds -> step}, {z, zp}, {h, hp}, MaxIterations -> 500]
]
We track the root:
solCurve1 = NestWhileList[
Check[predcorr1[#1, #2, #3], {#1, #2, #3/2}] & @@ # &, {{z0, h0},
{dz0, dh0}, 0.1}, #[[1, 2]] < 10 &] // DeleteDuplicates[#, (#1[[1]] == #1[[2]] &)] &;
Here is the bifurcation diagram:
ListPlot[Reverse /@ solCurve1[[All, 1]], AxesLabel -> {h, z}]
NSolve. That said, I'm not sure if your statement was in reference to this particular example or to a more difficult family of problems. $\endgroup$