ON THE NUMERICAL CONTINUATION OF ISOLAS OF EQUILIBRIA

We present a numerical strategy to compute one-parameter families of isolas of equilibrium solutions in ODEs. Isolas are solution branches closed in parameter space. Numerical 
continuation is required to compute one single isola since it contains at least one unstable segment. We show how to use pseudo-arclength predictor-corrector schemes in order to follow an entire isola in parameter space, as an individual object, by posing a suitable algebraic problem. We continue isolas of equilibria in a two-dimensional dynamical system, the so-called continuous stirred tank reactor model, and also in a three-dimensional model related to plasma physics. We then construct a toy model and follow a family of isolas past a fold and illustrate how to initiate the computation close to a formation center, using approximate ellipses in a model inspired by the Van der Pol equation. We also show how to introduce node adaptivity in the discretization of the isola, so as to concentrate nodes in region with higher curvature. We conclude by commenting on the extension of the proposed numerical strategy to the case of isolas of periodic orbits.


Introduction
Bifurcation theory for dynamical systems (see, e.g., [23]) is concerned with determining and classifying the long-term behaviour of evolution equations upon parameter variation. Continuous-time examples include ordinary differential equations, partial differential equations and delay differential equations. In this paper, we will consider problems involving smooth autonomous vector fields (1)ẋ = f (x; p, λ), x(t) ∈ R n , p, λ ∈ R.
The above system of equations describes the evolution in time of some physical or biological system depending upon a set of control parameters. Of particular interest are branches of attractors such as equilibria or limit cycles. An isola is generically defined as a closed curve in parameter space: points on an isola belong to the same class of invariant set, hence the curve is homotopic to a circle. We can have, for instance, isolas of equilibria or isolas of limit cycles. Once an isola is found by varying one of the bifurcation parameters, it is often interesting to study its persistence when a secondary parameter is changed. Keller ([20]) showed that Figure 1. Schematic representation of two families of isolas originating from two separate isola centers C 1 at λ = λ c1 and C 2 , at λ = λ c2 . Isolas z 1 (s) and z 2 (s) collide and give rise to a larger isola z 3 (s) at λ = λ c3 . We can continue z 1 for λ c1 < λ < λ c3 , z 2 for λ c2 < λ < λ c3 (two independent continuations spanning M 1 and M 2 respectively) and z 3 for λ < λ c3 (spanning M 3 ).
if the vector field f (x; p, λ) is smooth, then isolas occur in one-parameter families emanating from singular points, the isola (formation) centers (see also Figure 1).
Isolas of equilibria are well understood via singularity theory [16]; mechanisms through which an isola is born and annihilated have been classified systematically. Examples of isola destruction include mushrooms (when an isola collides with an unbounded branch) and broken pitchforks. These mechanisms have been observed, for example, in chemical kinetics, showing that isolas are important dynamical objects [15]. Isolas of periodic orbits are also dynamical objects to be expected in parametrized families of vector fields (see [33]). However, contrary to the equilibrium case, the underlying mathematical theory is not complete yet; in particular, the mechanisms behind the termination of such families of isolas are not fully understood.
In biological and physical applications the computation of a single isola is done via numerical continuation in the primary parameter p and there are several possible strategies to study its dependence upon the secondary parameter λ. In most cases, the exploration in λ is achieved by doing a series of p-continuations for slightly different values of λ (see [14]). Furthermore, since an isola contains at least two fold points, some authors continue these folds in two parameters to approximate the region of existence of isolas in the (p, λ)-space (see [18,14]).
Computing many isolas is tedious and requires an adequate (possibly systemdependent) strategy: it would be desirable to have a robust procedure to continue the isola automatically in any secondary parameter. The motivation behind the present paper is to propose a simple algebraic formulation that accomplishes this task and that can be easily integrated into existing continuation packages such as Auto [13], MatCont [11] and CoCo [8].
The main result of the present paper is that it is possible to compute branches of isolas of steady states using pseudo-arclength continuation, starting from a suitable initial guess. In this way, we can follow isolas towards both their birth (isola formation centers) and their death (mushrooms, broken pitchforks). In the former case, we can monitor the arclength of the isola in the proximity of the isola centers, and we can make use of previously developed techniques for the detection of the formation point [18,27]. In the latter case, the continuation will approach the branch point associated with the isola destruction (without reaching it), in a similar fashion to what happens when a branch of limit cycles moves towards a homoclinic connection. The main idea behind the proposed algorithm is to interpret isolas as level set curves (see also Figure 1). An isola is approximated by a closed polygon and the corresponding discretized problem is solved numerically via Newton method.
Numerical continuations of closed curves is usually done in connection with invariant circles of Poincaré maps for periodically-forced ODEs (non-autonomous systems similar to (1)). In that context, an invariant curve s → x(s) ∈ R n is computed as a fixed point of the Poincaré map Φ : R n → R n and various methods have been proposed to solve this problem (we refer to [31] for a more detailed introduction). A first algorithm, based on direct iterations of Φ and able to capture only attracting invariant curves, was proposed in [1]. Chan and Doedel, in [5], approximated the invariant curve with Fourier modes and applied collocation to the equation Φ(x(s)) = x(s + σ), fixing σ ∈ R via a suitable "anchor condition" (phase condition). In [22], the authors suggest a particular parametrisation of x(s), discretize it at some nodes s i and then impose Φ(x(s i )) ∈x(s), wherex(s) is a curve interpolating x(s) with interpolation points s i . In [31] and [32], invariant curves are approximated by a polygon through s i , and the points are carefully distributed, so that the polygon through Φ(x(s i )) converges to x(s).
More recently, a polygonal representation of closed curves was used to trace Arnol'd tongues [28]. As we shall see, isolas can be computed as zeroes of a suitable functional equation that can be discretized using Schilder & Peckham's approach.
We point out that the method proposed here is not the only way to explore attractors in higher-dimensional parameter spaces. In [7], the authors use numerical continuation to compute disjoint open branches of equilibria. Furthermore, a systematic approach to multi-parameter continuation was proposed in [17] and implemented in the package Multifario. Our strategy stands between a set of p-continuations for various values of λ ( [14]), and the full generality of the covering algorithm of Multifario. With respect to the former, we gain the possibility of spanning the (x, p, λ) space without writing scripts to switch from one continuation to the next, which becomes laborious when dealing with a large number of families of isolas; we also obtain, as a byproduct of our continuation, the arclength L of the isola, a quantity that varies considerably along the family of isolas, and that should be accounted for when designing the p-continuations mentioned above. With respect to the latter, we lose the ability to go past mushrooms or branch points within a single continuation.
The paper is organised as follows: in Section 2 we discuss how to set up the isola continuation problem; in Section 3 we show how to discretize the problem and we discuss our numerical implementation; in Section 4 and 5 we continue isolas in two applications, a continuous stirred tank reactor and a model from plasma physics, respectively; in Section 6 we show that our strategy can continue families of isolas that bend around a fold; in Section 7 we present an example in which isolas exist between two formation centers and comment on the computation of initial guesses; in Section 8, we show how to implement node adaptivity; in Section 9, we discuss possible extensions of the method.

Problem formulation for the continuation of isolas of equilibria
As stated in the introduction, we will focus on autonomous dynamical systems of the form expressed in Equation 1. The system might depend in general on more than two parameters. Without loss of generality, we consider here only two, namely p and λ, our primary and secondary parameter, respectively.
where the rectangular Jacobian matrix f z (z) ∈ R (n+2)×n has rank n for all z ∈ M. Furthermore, assume that M is homeomorphic to a cylinder.
Isolas of equilibria can be identified with closed paths on the manifold M. More precisely, we are interested in smooth closed curves s → z(s) = (y(s), λ) that lie in M for all s ∈ [0, 1], where we defined (x(s), p(s)) = y(s) ∈ R n+1 for convenience (see Figure 1). We can then formulate the computation of a single isola as follows: where the prime denotes differentiation with respect to s and · 2 is the Euclidean norm on R n+1 .
The computation of an isola can then be seen as an implicit Differential-Algebraic problem [2]. The first n equations constrain z(s) = (y(s), λ)) to lie in M, while the (n + 1)th equation states that, as we move on the manifold, the absolute value of the velocity remains constant. The scalar L is the arclength of the isola, as we have 1 0 y (s) 2 ds = L. We note that problem (2) has n + 2 unknowns (y(s), L), but only n + 1 boundary conditions. In fact, if y(s) is a solution to (2), so is y(s + q), for all q ∈ R. This arbitrariness can be eliminated by imposing an (n + 2)nd boundary condition: a similar problem is met in the numerical continuation of periodic orbits and it is usually solved by means of an integral phase condition ( [12,4,6]).
where y * is a reference solution (typically the solution of the previous continuation step) and ·, · denotes the standard inner product on R n+1 . Therefore, we propose the following algebraic problem for the computation of isolas.
Problem 2. (Isola computation) Assume that Hypothesis 1 holds. For fixed λ ∈ R, compute an isola as a solution (y, L) ∈ C 1 ([0, 1], R n+1 )×R + to the operator equation Remark 1. The problem (4) allows to continue isolas of steady states under the assumptions that M is sufficiently smooth. As mentioned in the introduction, isolas emanate from isola centers and are destroyed at branch points, therefore we expect that, as λ is varied, Hypothesis 1 is violated at critical points. An example is sketched in Figure 1. In general, we should expect to run multiple continuations and approach critical points from different directions.

Numerical discretization
Problem (4) can be discretized into a set of nonlinear algebraic equations and continued in λ using standard path-following routines as implemented, for instance, in AUTO [13]. To this end, we can choose a spline collocation method providing an approximate solution and a partition s 1 < s 2 < . . . < s m+1 of the interval [0, 1] which determines our collocation points.
We use m evenly-spaced nodes in [0, 1], that is, s i = (i − 1)h = (i − 1)/m and we approximate the isola via a closed polygon (see the schematic in panel (a) of Figure 2) with interpolation points given by In this setting, the discretized version of G is With the first mn components, we impose that each node is an equilibrium of the dynamical system (1); the following m components guarantee that nodes are uniformly spaced along the isola (with distances measured by the Euclidean norm). The perimeter l of the polygon is an unknown of the problem and the last equation of (5) is a discretized version of the integral phase condition (see also panel (b) in Figure 2). Since the isola computation can be seen as an interpolation problem on m nodes, convergence is guaranteed by the Weierstrass theorem [29], and we expect that l → L as m → ∞. Furthermore, the interpolation error is proportional to the second derivative of y(s), as we are using piecewise-linear interpolation (in Section 8 we show how to prescribe a node distribution that is denser where y (s) . . .
. . . is larger). We note that, together with the discrete isola (y 1 , . . . , y m ), we obtain simultaneously an approximation of all the isola folds. The discretization proposed here is similar to the one presented in [28, Section 2.2], where resonance surfaces are found by continuing closed curves in a secondary parameter. More precisely, we recover (5) by setting time derivatives to zero and removing boundary and phase conditions. We implemented pseudo-arclength isola continuations in Matlab, modifying EPCont (one of the ancestors of CoCo). In particular we amended EPCont so as to solve the nonlinear equations with fsolve, Matlab's routine for optimisation problems. In our implementation, we form explicitly the sparse Jacobian matrix associated with G and we use a trust-region reflective algorithm, setting the tolerance to 10 −6 .
In the remainder of the paper, we use our strategy in various examples with isolas of equilibria, which we follow from their birth to their annihilation. The computations shown in the next sections, where isolas are approximated with hundreds of nodes, run in about a minute on a standard laptop.

The continuous stirred tank reactor: isolas and mushrooms
The continuous stirred tank reactor (CSTR) is a well-studied planar model of chemical reactions (see [16,30]). The equations are We fix the following parameters: B = 8.0, β = 1.0, v c = 0.0. Initially, we set λ 0 = 1.1 and do a p-continuation to compute an initial isola and obtain an initial guess (y 0 , L 0 ) for the λ-dependent problem G(y, L; λ) = 0.
In Figure 3, we show a family of isolas obtained by continuation. The starting isola is marked in red, and we continued it, both for decreasing and increasing values of λ. As mentioned in the introduction, we expect the family of isolas to terminate in one of the following ways: p λ x Figure 3. Surface of isolas of the continuous stirred tank reactor in the (p, λ, x ) − space, with x = (u, v). We compute isolas with m = 500 evenly-spaced nodes.
(1) The arclength of the isola tends to zero as we approach an isola formation center. (2) The branch disappears at a branch point.
The family of isolas in Figure 3 features both scenarios: for decreasing λ, the isola shrinks down to a point, an isola formation center, that we located at approximately (p c , λ c ) = (0.23, 0.096). This value is obtained by inspecting the perimeter of the approximated isola, which can be monitored during our continuation (see Figure 4, panel b). For increasing λ, the isolas expand up to a point in which we observe a transition towards a mushroom, that is, an unbounded bifurcation curve containing a round-shaped part resembling an isola. This scenario is shown in Figure 4: in panel (a), we plot the family of isolas (in blue) as they approach the mushroom (in black), which has been computed with a standard p-continuation. As expected, the adopted algebraic problem can not go past the branch point, but we show that we can compute isolas arbitrarily close to it.
In general, we can continue isolas in any of the parameters λ, p, B, β, v c , of the CSTR system (6). More precisely, it is possible to start from an isola computed via the λ-continuation described above, and use it as the starting point for a continuation in any other parameter. In Figure 5 we compare a β-continuation (blue) and a λ-continuation (red) of the initial isola y 0 , plotted in the three-dimensional (p, u, v)-space. Similar results can be obtained in other continuation parameters.

A model from plasma physics: isolas and symmetry-breaking bifurcations
The second example that we used to test our continuation strategy is a threedimensional system modelling turbulent two-dimensional fluid motion in plasmas. The main differences with the previous example are the increased dimensionality of  the model, and the mechanism for the isola destruction: in this case, the isolas disappear through the break-up of two connected pitchfork bifurcations. The system is given by the following equations (for details about this model, see [3]): In our computations, we fix the parameters κ = 1.5, γ = 1, α = 2.4, β = 1, a = 0.3, b = 1. As usual, p and λ are our primary and secondary bifurcation parameters, but in this case the latter is related to the Z 2 -symmetry of the model, as system (7) is left invariant by the transformation (w, λ) → (−w, −λ).
In Figure 6 we show branches of equilibria in the parameter p, for λ = 0 (black curves). The trivial branch corresponds to the zero solution and gives rises to two asymmetric branches, that emerge and disappear via a subcritical (PF 1 ) and a supercritical (PF 2 ) pitchfork bifurcation. We expect that this bifurcation diagram will not be robust to perturbations in λ: by repeating the p-continuation for λ = 0.02 (blue curves in Figure 6), we find an isola and a coexisting open branch, as each pitchfork bifurcation of the unperturbed system gives rise to two fold bifurcations (not labelled in the figure).
Similarly to what was done in the CSTR system, we interpolated the isola at λ = 0.02 on m = 300 nodes and continued it for both increasing and decreasing value of λ. As λ → 0 + , we approach the unperturbed configuration, and the topleft and top-right fold points converge towards PF 1 and PF 2 respectively. We then continue from λ = 0.02 for increasing λ and find a family of isolas terminating at an isola formation center at about (p c , λ c ) = (15.76, 2.43). The complete family of β-continuation λ-continuation p u v Figure 5. Two parametrised families of isolas of the CSTR system (6) are represented in the state/parameter (u, v, p)-space. A β-dependent branch is plotted in blue and a λ-dependent branch, starting from the same y 0 , is plotted in red. In this computation we set m = 500.
isolas obtained for λ > 0 is shown in Figure 7 as a blue surface, labelled I + , in the (w, p, λ)-space.
By symmetry, we infer the existence of a second family of isolas, I − , symmetric to I + with respect to the axis λ = 0. The new family is also found by continuation, and it is plotted in red in Figure 7. In panel (b), we also show a few unbounded branches existing for λ = 0 and corresponding to deformations of open structures found in Figure 6 for λ = 0.02.

Continuing a family of isolas past a fold
One of the main advantages of pseudo-arclength continuation is that it allows to compute solution branches of nonlinear equations past a fold point. In the context of stationary and periodic solutions of parametrised vector fields, this provides information about stable and unstable families of such attractors, and fold points of the solution branch correspond to saddle-node bifurcations of the underlying vector field. Isolas are not dynamical objects so there is no sensible concept of stability attached to them. They are formed by a parametrised family of attractors, each of them having a specific stability -note that this stability can change multiple times along an isola. Nevertheless, as any solution branch to a set of nonlinear equations, . One-parameter bifurcation diagram of system (7) in the parameter p. For λ = 0 (black curve) we find two pitchfork bifurcations (black dots, labelled PF 1 and PF 2 respectively) and for λ = 0.02 we find an isola (blue curve).  Figure 6 (blue surface, labelled I + ). By symmetry, we infer the existence of another family (red surface, labelled I − ) occurring for negative values of λ; panel (b) also shows unbounded branches that coexist with the isolas for λ = 0. We approximate each isola with m = 300 nodes.
a family of isola can have folds, and we can compute it using our path-following strategy.
To illustrate this point, we construct a toy model of a one-dimensional dynamical system with two parameters, such that its two-parameter family of equilibria lies on a folded surface that looks locally like a Mexican hat (see Figure 8 (a)). The toy model is given by where p and λ are the primary and the secondary parameter, respectively. We fix: a = 0, b = 1. In this very simple example, we know directly from equation (8) that the family of equilibria is the graph {(x, p, λ) ∈ R 3 | λ = Γ(x, p)} of the function Level sets of Γ correspond to a p-dependent family of equilibria for a given fixed value of λ. From equation (8), it is clear that we have an isola or two coexisting isolas for values of λ away from the critical points of Γ, that is, for When λ = ab, that is, for x = p = 0, the function Γ has a non-degenerate critical point and the unique equilibrium is the origin. This corresponds to the isola formation center (red dot in Figure 8 (b)) where Γ x = Γ p = 0 and the Hessian is non-degenerate; a family of λ-dependent isolas emerges, initially with arbitrary small diameter. When λ = −(a − b) 2 /4, the function Γ has a ring of degenerate critical points (non-Morse [25]) given by x 2 + p 2 = (a + b)/2. This corresponds to the fold of isolas (bold red in Figure 8), where two family of isolas coalesce and disappear when decreasing λ further, for our choice of a and b (Γ x = Γ p = 0 and the Hessian is degenerate along this special isola). The result of our continuation of isolas is presented in Figure 8. Panel (a) shows a three-dimensional view of the family of isolas in parameter/variable space (x, p, λ), close to the fold (highlighted in red), for m = 100; panel (b) shows a twodimensional projection onto the parameter plane (λ, p). We initiated our computation with a small isola obtained with traditional continuation routines for λ = −0.2. Subsequently we continued the initial isola using the proposed formulation in both directions. For decreasing values of λ, the inner family of isolas were continued p x λ Figure 9. Family of isolas of system (9)-(10) for |λ| < 1. The starting isola is shown in bold red, the continuation is done for decreasing lambda down to the second isola center at λ = −1 (black dot). Each isola is approximated with m = 300 nodes. towards the formation center at the origin; in the opposite direction, we continued the inner family past the fold at λ = −(a − b) 2 /4, towards the outer family, which extends up to infinity.

Computing initial guesses in the proximity of isola centers
So far, we started our continuations from a previously-computed isola. Even though an initial guess might not be available, one can always search for a formation center and find a nearby isola with a small perimeter. The detection of isola formation centers can be done via solving an appropriate extended system during a two-parameter continuation (see [20,27] and [21] for a detailed discussion on the implementation details, including continuation of isola centers in a third parameter).
We now illustrate the finding an initial guess when an isola formation center has been detected. To this end, we follow the analysis developed in [9], where the authors give equations for ellipses approximating isolas close to the bifurcation centers. The following example model is designed so that it possesses a family of isolas existing in a compact interval of the secondary parameter λ, each end corresponding to a formation center. The resulting isola family has then the shape of an ellipsoid (see Figure 9). The system is two-dimensional and its equations arė To construct system (9)-(10), we took the fast equation of a Van der Pol type model and shifted it by a constant p to obtain the second equation. Then one can easily eliminate y from the equilibrium condition and rewrite it as the quadratic equation the solutions of which are given by From equation (12), we can deduce that two real solutions (or a double solution) exist provided the condition p 2 ≤ 12(1−λ 2 ) is satisfied. This means that for |λ| < 1, equation (12) defines an ellipse in the (p, x)-plane. Furthermore, this ellipse shrinks to a single point at λ = ±1 and there is no real solutions for |λ| > 1. Therefore, we can conclude that system (9)-(10) possesses a p-dependent family of isolas of equilibria that is created and destroyed via isola centers at λ = ±1.
We now consider the initiation of an isola continuation close to the formation center at λ = 1. Following [9], we take an ε-perturbation λ = 1−ε/2 (for 0 < ε 1) and we compute approximations to the ε-dependent family of isolas of system (9)- (10). To this end, we expand the left-hand side of (11) about the isola center. We then truncate the resulting expansion to second order in ε, which yields the following equation for the ellipses Equation (13) corresponds exactly to equation (11) for λ = 1 − ε/2, hence for the case considered here the initial guess is also an exact solution. Note that when the equilibrium equation is vectorial, the perturbation method used above is more elaborate but works equally well and can be implemented in any continuation package; we refer the reader to Section 2 of [9] for more details.
The result of the isola continuation in system (9)-(10) is presented in Figure 9; the computation is initiated with an ellipse given by equation (13) for ε = 0.01, that is, λ = 0.995. We show the entire family of isolas starting with the bold red isola at λ = 0.995, very close to the right isola center (λ = 1), down to the other isola center at λ = −1. We compute 60 isolas to render this bubble-shaped family, with m = 300 nodes on each isola.

Node adaptivity
The problem (5) offers the possibility to implement adaptivity in the node distribution. Since we have approximated the isola with piecewise-linear polynomials, the interpolation error is proportional to the curvature of the isola. As we have seen in the examples treated in Sections 4 and 5, isolas may deform along the computed family and develop regions of high curvature. A possibility to overcome this problem is to use a large number of nodes: while this is computationally affordable in the examples we considered so far, it is desirable to introduce node adaptivity for more demanding problems. This can be achieved by discretizing the continuous operator G in a different way: instead of distributing the nodes uniformly along the isola (with respect to the Euclidean distance), we distribute uniformly the interpolation error. This amounts to replacing equations  Figure 10. Illustration of the node adaptivity. We distribute the interpolation error evenly along the isola. The curves are isolas of the toy model (14) for λ ∈ [0.1, 2.8], discretized using m = 120 nodes.
where e(y i−1 , y i , y i+1 ) is the local interpolation error and σ is a new unknown of the problem. In the case of piecewise-linear interpolation, the error can be approximated by therefore the nodes adjust so as the interpolation error remains constant at the fixed value σ [26]. We point out that more sophisticated adaptations are also possible. In [28], the authors propose a hybrid between uniform node distribution and uniform errors distribution. We illustrate the effect of node adaptivity with the toy model for which we continue isolas for λ ∈ [0.1, 2.8] using m = 120 nodes. In Figure 10, we show the projection of the family of isolas onto the (p, x)-space. As is clearly illustrated in this figure, the nodes change their position during the continuation, accumulating in regions with higher curvature.

Conclusions
We presented a numerical continuation strategy to compute one-parameter families of isolas of steady states. The main idea behind our method is to interpret isolas as individual objects: each isola is approximated by a polygon, and for fixed values of the secondary parameter we solve an algebraic problem whose unknowns are the nodes of the polygon and its perimeter. We fix the phase of the isola with an additional equation, similar to that used in the continuation of periodic orbits. The nodes distribution can be uniform (equally-spaced nodes, with distance measured by the Euclidean norm) or adaptive (distributing evenly the interpolation error).
We applied our method to continue isolas of equilibria in several examples, namely the continuous stirred tank reactor, a model from plasma physics, a toy model displaying a fold of isolas, and a fourth model inspired by the Van der Pol equation and displaying a family of isolas bounded by two formation centers. In all the test cases, we show that it is straightforward to set up a continuation problem, compute an initial isola, interpolate it on a large number of nodes, and continue it in parameter space.
The continuation proposed here is based on the fundamental assumption that isolas of equilibria lie on a sufficiently smooth manifold. Typically, families of isola terminate either at an isola formation center, or at a branch point. They correspond to cases where our initial assumption is violated and they need a special treatment. Isola centers can be detected (and continued) accurately by using the algorithms proposed in [27,30]. Furthermore, if we need to start the continuation from an isola center and there is no other indication of a nearby isola, it is possible to find a good initial guess by using the method proposed by [9].
On the other hand, the proposed method does not allow to follow isolas past mushrooms or branch points with a single continuation run. We remark that this should not be seen as a limitation, since the family of isolas terminates at such critical points (with the exception of the case outlined in Figure 1), However, it is always possible to start a new continuation past the criticality, in the same spirit of what is currently done in standard continuation packages when similar circumstances occur. Passing from M 1 to M 3 in Figure 1, for instance, is similar to branch-switching near a period-doubling bifurcation (we use the same algebraic formulation and start from an appropriate initial guess past the criticality). Approaching a mushroom (see Figure 4) is similar to continuing periodic orbits in the proximity of a homoclinic connections (for which one has to resort to a different boundary-value problem).
Since our formulation can be implemented for generic smooth vector fields, this strategy can be included in any existing continuation package (this extends to packages that can continue equilibria of delay-differential equations, for which the isola continuation reduces to a problem of the type (5)). The computational cost of the proposed strategy amounts to discretize an additional spatial direction s, and it is feasible in the case of isolas of equilibria, as our numerical experiments suggest.
Finally, the formulation can be extended to continue isolas of periodic orbits, in which each node on the isola is determined by a parameter value p i , a periodic orbit x i (t), t ∈ [0, 1], and the relative period T i . For each node, we can replace the equilibrium condition with a boundary-value problem of the form The corresponding condition for the node distribution can be written as d(y i , y i+1 ) = y i − y i+1 , where the Euclidean norm in (5) is replaced by a suitable norm on C 1 ([0, 1], R n ) × R [28]. In this case, an integral phase condition for the isola can also be defined in order to close the problem. However, owing to the large number of unknowns involved in this case, it would be necessary to implement an adaptive algorithm, where the adaptivity is to be intended both along the isola (accumulating the nodes where the curve becomes sharper) and in the discretization of the periodic orbits (each node having a separate set of collocation points).