Spiraling bifurcation diagrams in superlinear indefinite problems

This paper computes and discusses a series of intricate global bifurcation 
diagrams for a class of one-dimensional superlinear indefinite 
boundary value problems arising in population dynamics, 
under non-homogeneous boundary conditions, measured by $M>0$; the 
main bifurcation parameter being the amplitude $b$ of the superlinear terms.


1.
Introduction. This paper computes the bifurcation diagrams of positive solutions for a series of symmetric and non-symmetric one-dimensional boundary value problems −u = λu + a b (x)u p in (0, 1), where M > 0, p > 1 and λ < 0 are constants and a b (x) is a weight function that changes sign in (0, 1), of the form where α ∈ (0, 0.5), b > 0 and c > 0 are constants, and χ ∈ C[α, 1−α] is a continuous function such that χ(x) > 0 for all x ∈ [α, 1 − α]. Superlinear indefinite problems of this type, mostly when M = 0, have attracted a great deal of attention during the last two decades, as it is apparent by having a look at Berestycki, Capuzzo-Dolcetta and Nirenberg [4,5], Alama and Tarantello [1], Amann and López-Gómez [3], Gómez-Reñasco and López-Gómez [22,23], Mawhin, Papini and Zanolin [42], López-Gómez [28][29][30], García-Melián [18], and the list of references therein. When M = ∞, these solutions are referred to as the positive large solutions of In the sublinear case when b ≤ 0 and M < ∞, it is well known that Problem (1) possesses at most one solution, which is a global attractor for    ∂ t u − ∂ xx u = λu + a b (x)u p in (0, 1), t > 0, u(0, t) = u(1, t) = M, t > 0, u(·, 0) = u 0 > 0, in (0, 1), if it exists (see López-Gómez [31]). However, when b > 0, this situation changes dramatically. Indeed, for the simplest choice it has been recently discovered that the number of positive solutions of (1) increases arbitrarily as λ < 0 decreases approximating −∞ (see [41]). Subsequently, a solution u(x) of Problem (1) is said to be of type n ≥ 0 if it has n strict critical points (either local maxima or minima) in the interval (α, 1 − α). More precisely, u is said to be of type (n, a) if u is of type n and asymmetric around x = 0.5, while it is said that it is of type (n, s) if it is of type n and symmetric. The main available multiplicity result for Problem (1) with the choice (4) is Theorem 3.6 of López-Gómez, Tellini and Zanolin [41], which can be summarized as follows Theorem 1.1. Fix c > 0, p > 1 and M > 0. Then, there exists a decreasing unbounded sequence {λ j } j≥1 , with λ 1 = 0, such that, for every j ≥ 1 and λ ∈ (λ j+1 , λ j ), there exists b * = b * (λ) such that Problem (1) for the choice (4) with b = b * admits, at least, 2j solutions; one of them of type 0, another of type 1, and two of type n, for each 2 ≤ n ≤ j, if j ≥ 2. Moreover, all the solutions of type n = 2k, k ≥ 1, are asymmetric in [0, 1], while the remaining ones are symmetric.
Actually, for every j ≥ 1, one has that λ j+1 ∈ τ −1 j (1 − 2α), where the τ j 's are the Poincaré maps defined in [41,Section 3]. Furthermore, perturbing globally from b = b * one can obtain all the possible global bifurcation diagrams of positive solutions of Problem (1) for the choice (4). Essentially, a series of loops bifurcate sequentially from a certain primary branch that spirals around the zero-type solution as λ < 0 decreases; the number of loops increasing arbitrarily as λ decreases (see Figures 1,2,4,5,6,7 and 9). More precisely, after every half rotation of the primary curve around the zero-type solution, a new loop emerges bifurcating from it. As a consequence, the following result holds. x ∈ [0, α), ∞, x ∈ [α, 1 − α], u s (1 − x), x ∈ (1 − α, 1], where u s stands for the unique solution of the singular problem −u = λu − cu p in (0, α), u(0) = M, u(α) = ∞, plus, possibly, finitely many secondary loops, emerging each one at two symmetrybreaking bifurcations points from the primary curve, which are filled in by asymmetric solutions. Actually, as λ < 0 decreases, the primary curve spirals over itself and the number of emanating loops connected to it increases arbitrarily. More precisely, with the notations of Theorem 1.1, if λ j+1 < λ < λ j for some j ≥ 1, then the minimal bifurcation diagram possesses • at least j turning points on the primary curve and n loops bifurcating from it (with the possibility of having an additional secondary loop) if j = 2n + 1 with n ≥ 0. • at least j −1 turning points on the principal curve and n loops bifurcating from it (with the possibility of having 2 additional turning points on the principal curve) if j = 2n with n ≥ 1.
Rather astonishingly, though the simplest prototype model can have an arbitrarily large number of positive solutions, the following uniqueness result, going back to López-Gómez, Molina-Meyer and Tellini [39] holds. It is remarkable that it is valid independently of the choice of χ in any spatial dimension. In population dynamics, the solution u(x, t) of (3) models the density of the species at the location x ∈ (0, 1) at time t > 0, λ < 0 measures the net death rate of the species, and u 0 > 0 is the initial population density. In these models, λ < 0 when pesticides are used in high concentrations, or a certain patch of the environment is polluted by introducing chemicals, waste products, or poisonous substances. In spatial ecology, the function bχ(x) measures the interspecific facilitative effects of the species u in the patch (α, 1 − α). In these generalized logistic prototypes, the individuals of the species u compete for the natural resources in the region where a b < 0, while they cooperate in the patches where a b > 0. Although there is a strong experimental evidence of interspecific competition (see, e.g., Shoener [46] and Connell [13]) and positive interactions are well documented among organisms from different kingdoms (see, e.g., Hutchinson [24], Wulff [47] and Saffo [45]), finding positive interactions between similar organisms seems to be a huge task in empirical studies, as they do not arise alone but in combination with competition. However, according to the abiotic stress hypothesis of Bertness and Callaway [6], the importance of positive interactions in plant communities increases with abiotic stress or consumer pressure. Several empirical studies support the validity of the abiotic stress hypothesis, for example a substantial amount of positive interactions in plant communities have been isolated in harsh environmental conditions (see, e.g., Callaway, Walker [10] and Pugnaire [44]). According to these monographs, (3) is a reasonable model for studying the effects of combined facilitation and competition in polluted habitat patches. Theorems 1.1 and 1.2 reveal that, under facilitative effects in competitive environments, the harsher the environmental conditions, the richer the dynamics of the species. Being this a rather paradoxical result, the fact that the number of solutions of (1) grows arbitrarily, as an effect of interspecific facilitative effects, in spite of the growth of the degree of hostility in the habitat, might be a significant feature in the theory of ecosystems explaining the biodiversity of the Earth.
To analyze the dynamics of (3), whose steady states are the positive solutions of (1), it is imperative to design efficient numerical algorithms to compute the positive solutions of (1), as well as the dimensions of their unstable manifolds, because they cannot be ascertained, nor even estimated, with the available abstract analytical results. This paper computes all the underlying global bifurcation diagrams of (1) for several choices of χ(x), by using some global path-following solvers in b; λ playing the role of a secondary bifurcation parameter.
By some classical results in elliptic regularity theory (see, e.g., Gilbarg and Trudinger [19]), any weak solution u of (1) must satisfy Moreover, thanks to the strong maximum principle (see, e.g., [34]), any positive solution of (1) must satisfy u(x) > 0 for all x ∈ [0, 1]. The solution with 0 strict critical points, referred to as the trivial solution of (1), will play an important role in explaining and understanding the complexity of the global bifurcation diagrams of (1) with choice (4), as it plays a similar role to an organizing center in singularity theory (see Golubitsky and Shaeffer [21]). To construct it, let u (x) be the unique solution of −u = λu − cu p in (0, α), u(0) = M, u (α) = 0, and set m 0 := u (α). Then, the constant function m 0 solves , and, for such value of b, the trivial solution u t is defined through Naturally, u t is symmetric around x = 0.5. Basically, as −λ > 0 increases, a piece of the primary curve rotates counterclockwise around the trivial solution u t and, almost after every half rotation, an additional closed loop emanates from this curve. The loop consists of solutions of asymmetric type and it persists for all further values of −λ (see Figures 1,2,4,5,6,7 and 9).
Besides confirming and illuminating all the previous available analytical results, our numerics reveal some new interesting features not previously observed in the literature. Among them, the astonishing phenomenon that lim λ↓−∞ b c (λ) = ∞ (see Theorem 1.3), whereas, simultaneously, lim λ↓−∞ u = 0 uniformly in [0, α]∪[1−α, 1], which entails that all the turning points along the primary curve are extremely closed. Figure 9 illustrates very well this phenomenology. It shows a plot of the global bifurcation diagram of positive solutions of (1) for λ = −2000. Although b c (λ) ∼ 1.2463985×10 7 , most of the positive solutions, independently of their types, satisfy u(α) ∼ 10 −4 . Moreover, the first closed loop bifurcates from the primary branch at b ∼ 1.2463984 × 10 7 , which is extremely close to the turning point b c (λ). This strong squashing effect of the diagrams as λ ↓ −∞ makes the adaptation of the available continuation codes to carry out the numerical computations of this paper a very difficult task.
Another astonishing phenomenon is the fact that the component of the set of solutions of (1) emanating from the unique positive solution of (1) at b = 0 for the special choice (4), splits into an arbitrarily large number of components as χ perturbs from 1 in an asymmetric way, because each of the loops of asymmetric solutions emanating from the primary branch generates an isola (see Figures 14,15 and 16). Consequently, the number of components might change dramatically when the weight function a b loses its symmetry, although the numerics show that similar multiplicity results to the one of Theorem 1.1 should be true for any bounded weight function a b (x), as soon as it changes sign.
As a consequence from our numerics, we were also able to determine the dimensions of the unstable manifolds of all positive solutions of (1). It turns out that in the special case χ = 1 they change by 1 every time a turning or a bifurcation point is crossed, increasing until we reach the trivial state, and then decreasing.
From the point of view of the applications in population dynamics, our results establish that, under facilitative effects in competitive media, the harsher the environmental conditions, the richer the dynamics of the species.
Besides the Introduction, this paper consists of three sections. Section 2 describes our numerical experiments for the special case χ = 1, Section 3 analyzes the changes in the structure of the bifurcation diagrams, both for symmetric and asymmetric weights, as χ perturbs from 1. Finally, Section 4 discusses the numerical codes used to carry out the numerics of this paper, focusing the attention on the main computational troubles provoked by the strong squashing effects of the diagrams. Therefore, it has an intrinsic interest from the point of view of sharpening the efficiency of the available path-following algorithms.

2.
A series of significant bifurcation diagrams for χ = 1. Throughout all this section we will assume (4) and we will fix p = 2, M = 100, α = 0.3, and c = 1, whereas b and λ will be regarded as the primary and secondary parameters of the problem, respectively. Since the value of c is substantially smaller than M , the associated solutions of (1) can be though of as approximations to the large solutions of (2).
Precisely, we will give to the parameter λ a series of significant values ranging in the interval (−∞, 0) and, for each of these values, b will be regarded as the main bifurcation parameter to compute the corresponding bifurcation diagrams of (1). Our main goal is to analyze their complexity as −λ > 0 increases. Figure 1 shows the plot of the global bifurcation diagram of (1) for λ = −5, on the left, as well as the profiles of a series of solutions along it, on the right.
As in all subsequent global bifurcation diagrams, in the left picture we are representing the value of b, on the horizontal axis, versus the value of u(α), on the vertical one. The bifurcation diagram consists of a single primary curve emanating from the unique classical positive solution of −u = λu + a 0 (x)u p in (0, 1), at b = 0, whose existence and uniqueness was established in López-Gómez [31], and it continues towards the right up to reach the critical value b = 0.3401, where it goes backwards exhibiting a subcritical turning point. Once passed the turning point, the solutions on the upper half-branch can be continued for every 0 < b < 0.3401. All our numerical experiments have revealed that along the upper half-branch the solutions of (1) blow up in [α, 1 − α] as b ↓ 0, while in [0, α) they approximate the unique solution u s of the singular problem (6). This means that, as b ↓ 0, these solutions approximate the metasolution m(x) defined by (5), as conjectured by López-Gómez, Tellini and Zanolin [41]. Therefore, the global bifurcation diagram establishes an homotopy between the unique classical solution of the sublinear problem (7) and the unique metasolution of (2) supported in (0, α) ∪ (1 − α, 1), defined in (5) (cf. López-Gómez [31,32]). Consequently, for every b ∈ (0, 0.3401), (1) admits, at least, two solutions. Moreover, the solutions along the lower half-branch are linearly asymptotically stable, while the solutions along the upper half-branch are unstable with one-dimensional unstable manifold.
The bottom plot of the right picture in Figure 1 shows the unique solution of (7), the next one shows the trivial solution u t , whose exact position in the global bifurcation diagram has been marked with a thick point, the third one shows the solution of (1) at the turning point and the fourth one shows the solution on the upper half-branch for b = 0.2097. The trivial solution arose at b = b * = 0.1811 on the lower half-branch. All the solutions on the lower half-branch to the left of b * are convex in the central interval (0.3, 0.7), and stay below m 0 , while the remaining solutions of the diagram are concave and stay above m 0 . A remarkable feature is that the solutions along the global bifurcation diagram are pointwise increasing as we run over it starting at the unique solution of (7). Therefore, they decrease with b along the upper half-branch, while they increase with b along the lower one; a genuine superlinear behavior. The solutions along the upper half-branch have been computed up to b = 2.5 × 10 −7 , where u(α) ∼ 8 × 10 5 and u(0.5) ∼ 3 × 10 8 .
Similarly, the left picture of Figure 2 shows the global bifurcation diagram for λ = −70 and the right one a series of solution plots along it. The unique significant qualitative difference with respect to the previous case is that for λ = −70 the trivial solution u t arises on the upper half-branch of the diagram, instead of on the lower one. Actually, as λ decreased from −5 up to reach −70, the trivial solution u t moved along the lower half-branch and crossed the turning point at some intermediate value of λ up to reach the position marked by the thick point in the global bifurcation diagram of Figure 2.
The plots of the right picture of Figure 2 show the solutions of (1) on the lower half-branch for b = 0, b = 7.6585 (the turning point), and on the upper half-branch for b * = 6.4367 and b = 3.4651, 1.2513 and 0.8720, from the bottom to the top. As in the previous and in all subsequent cases, the solutions on the upper half-branch do approximate the metasolution m(x) defined in (5), as b ↓ 0. As before, all the solutions filling in the lower half-branch are linearly asymptotically stable, while the solutions along the upper half-branch are unstable with one-dimensional unstable manifold.  Figure 3 shows a magnification of the (subcritical) turning points exhibited by the two previous global bifurcation diagrams. Having a careful look at the scales at which they have been represented, it becomes apparent that the turning point for λ = −5 is more open than the corresponding one for λ = −70, which goes backwards at a much faster rate. Such behavior has been observed for all values of λ for which the global bifurcation diagrams of (1) have been computed. The more negative the parameter λ, the faster the rate at which the turning point is crossed. This causes one of the main technical difficulties we had to overcome in order to compute all the turning points for λ ≤ −800, where such a phenomenon is extremely emphasized, as it will be described in Section 4. In the previous cases, all the solutions computed along the bifurcation diagrams are of type (1, s), except u t , which is of type (0, s). Actually, the solutions along the curve of the bifurcation diagrams exhibit a minimum before reaching u t and a maximum after passing it. A similar phenomenon, consisting in interchanging minima with maxima, occurs every time we cross u t on the primary curve for all λ < 0.
Further, as λ decreased from λ = −70 up to reach the value λ = −140.868 a solution loop emanated from the upper half-branch of the primary curve. When the loop bifurcated from the primary curve, all solutions on it were of type (1, a). Figure  4 shows a series of emerging loops for four decreasing values of λ. They have been plotted with a red dashed line. The first plot of Figure 4 shows a significant piece of the upper half-branch of the primary curve containing one of these loops, where we have also marked the trivial solution u t through a thick point. Naturally, since the remaining solutions along the global bifurcation diagram are of type (1, s), the loop emanated as a consequence of a symmetry breaking bifurcation phenomenon as λ decreases. The bifurcation points of these loops are of pitchfork type. As λ further decreased from λ = −140.868, the closed loop of asymmetric solutions gradually enlarged up to reach the trivial solution u t at some critical value of the parameter λ.     [41]. All the solutions along these closed loops are unstable, with one-dimensional unstable manifold, while all the solutions on the portion of the primary branch encircled by the loop have two-dimensional unstable manifolds. This fact is in complete agreement with the local bifurcation theorem and the exchange stability principle of Crandall and Rabinowitz [14,15]. These numerical results reveal moreover that the minimum of the Poincaré map θ 2 introduced in Section 4 of [41] does not vary monotonically with b, because the emerging loop in Figure 4 left outside the trivial solution u t when it bifurcated from the primary curve.
As λ decreased from λ = −150, there was a critical value of λ where a bifurcation to solutions of type (3, s) from u t occurred, in such a way that the solutions on the primary curve around u t were of type (3, s) before changing type to (1, s). Figure 6 shows the global bifurcation diagram of (1) for λ = −300 on the left and, on the right, the plots of the trivial solution and two solutions of type (3, s) along the primary branch. They were computed for the values b = 227.6725 and b = 353.6714. The global bifurcation diagram of Figure 6 consists of two curves: the continuous blue line, which is the primary branch connecting the unique solution of (7) with the metasolution (5), and the dashed red line, which is the closed solution loop emanating from the primary curve at the values of the parameter b = 12.8294 and b = 526.4099, which are the turning points of the bifurcated closed loop. The subcritical turning point of the primary curve occurs at b = 527.4319. All solutions on the primary curve between the two squares are of type (3, s), while the remaining ones are of type (1, s). The changes between these two types occur at the values b = 149.8536 and b = 409.4183. Observe that all the changes of type on the primary branch occur at the level of u t , and this is a general feature, as proven in [41]. The type of the solutions along the closed loop follows the general A non-expected and rather striking feature is that the subcritical bifurcation point of the loop from the primary branch is very close to the turning point of the primary branch, as shown in the central plot of Figure 6. Although in this case the proximity of the two curves does not give any special numerical difficulty, provided that the step of the path following solver is sufficiently small, the proximity of the bifurcation point and the subcritical turning point is a very subtle question that will be discussed in Section 4.
As λ decreased from λ = −300 up to reach the value λ = −750, whose associated global bifurcation diagram has been plotted in the first picture of Figure 7, the arc of curve of the solutions of type (3, s) along the primary branch rotated counterclockwise around u t originating two additional turning points along it. Naturally, one of them supercritical and the other one subcritical. In this case the trivial solution arises at b * = 1.4732 × For a certain value of λ < −800 a bifurcation to solutions of type (5, s) from u t occurred. Then, the arc of curve of the primary branch encircled by the second loop rotated counterclockwise around u t originating two additional turning points on the (5, s)-arc of curve surrounding u t . Once these turning points had arisen, they separated gradually as λ decreased. The last plot of Figure 7 shows the global bifurcation diagram computed for λ = −1300. Topologically, this bifurcation diagram is equivalent to the previous one, but it shows two additional turning points on the primary curve, consisting of solutions of type (5, s), which became of type (3, s) before abandoning the interior of the second loop.
The last plot of Figure 8 shows the graphs of two solutions of type (5, s) for λ = −1300 with b = 2.7238 × 10 5 , the one with three minima, and b = 3.3598 × 10 5 , the one with three maxima. Naturally, to switch from each of them to the other, one must cross the trivial solution u t along the primary curve.
Although the standard path following routines were enough to compute all the previous global bifurcation diagrams until λ = −760.3, they started to give a number of computational troubles for λ ≤ −800. In Section 4 we will discuss how we were able to overcome them. Figure 9 shows the global bifurcation diagram of (1) for λ = −2000. Now, the counterclockwise rotation around u t exhibited by the last diagram of Figure 7 has been tremendously magnified, and a new closed loop, the third one, emanated from the solutions of type (5, s) along the primary curve. In Figure 9 such loop has been represented through a dotted purple line. As usual, the small squares on it are marking the points where the solutions of type (5, a) changed their type to (6, a).
To explain the global bifurcation diagram of Figure 9, we will run over the entire primary curve to describe its main features. We started by computing the unique solution of (7). Then, we continued the curve for increasing values of b and computed the branch of minimal solutions, all of them linearly stable, until we met the first (subcritical) turning point, arisen at b = 1.2463985 × 10 7 , where the solutions became unstable with one-dimensional unstable manifold until we reached the bifurcation point of the first loop, at b = 1.2463984 × 10 7 . Then, the solutions along the primary curve became unstable with two-dimensional unstable manifold.
Up to here, all the computed solutions were of type (1, s). Next, we continued the primary curve until b = 1.0292 × 10 7 , where the type of the solutions changed An extremely remarkable feature is that the dimension of the unstable manifolds of all solutions along the primary curve increased by one at each bifurcation or turning point we crossed, until we reached the interior of the last emanated loop, where these dimensions began to decrease, according to the same rule, until they became one-dimensional again.
As far as the behavior of the n-th loop is concerned, either it fully consists of solutions of type (2n − 1, a), as it occurs soon after its bifurcation from the primary branch, or it consists of solutions of type (2n−1, a) near the bifurcation points from the primary curve together with two central arcs of solutions of type (2n, a), for each n = 1, 2, 3, . . . . In any circumstances, the dimension of the unstable manifold of the solutions on the n-th loop is 2n − 1.
Note that b = 1.2463985 × 10 7 , the value where the first turning point along the primary curve occurs, is extraordinarily high, not only with respect to the value c = 1, but also with respect to u(α) = 2.1950 × 10 −4 , whereas u(0) = u(1) = 100. As the quotient b/u(α) is of order 10 11 , we are definitely close to the limits of the machine precision. Consequently, our numerical code seems extremely efficient. Figure 10 shows It should be noted that, for every b in between the two bifurcation points of the third loop from the primary curve, Problem (1) possesses at least 12 solutions: 6 among them symmetric and the remaining 6 asymmetric.
For smaller values of λ, we were not able to automatize the continuation codes, since we were outside the machine precision range, as we will show in Section 4. So, we stop our discussion here.
3. Some significant bifurcation diagrams for more general χ's. This section provides some of the most striking results of our very last numerical experiments for a series of perturbations from χ = 1. Precisely we will consider weights a b corresponding to two types of functions χ(x) and whose plots have been represented in Figure 11. Figure 11. Plots of a b for χ = χ 3 (A) and χ = χ 3,2 (B).
We will first consider the symmetric functions with n ∈ 2N + 1, 0 < α < α 1 < 0.5 and 0 < µ < 1, so that χ n > 0 in (α, 1 − α). Note that n is the number of oscillations of χ(x) around 1 in [α 1 , 1 − α 1 ] (see Figure  11 (A)) and it has to be odd in order to have a symmetric function. Moreover, χ 0 = 1. Later, we will consider the asymmetric perturbations with n, α, α 1 and µ as before and ε > 0. Observe that the main effect of ε is to modulate the amplitude of the last hump of χ n and that χ n,1 = χ n . A paradigmatic plot of these asymmetric weights has been already represented in Figure 11 (B). In all the subsequent numerical experiments we took n = 3, α 1 = 0.45 and µ = 0.1, while the other parameters were left unchanged with respect to Section 2, in order to compare the corresponding diagrams with those already discussed there in. Moreover, like in Section 2, we are always representing the values of b, on the horizontal axis of the diagrams, versus the values of u(α), which are plotted on the vertical one. Figure 12 shows, in the left column, the plots of four bifurcation diagrams computed for different values of λ in the case χ = 1, and, in the right column, the plots of the diagrams for the same values of λ as in the left, but for the choice χ 3 , instead of χ = 1. In this case, as it can be easily realized by simply having a glance at the plots of Figure 12, for the same value of λ, there is no significant qualitative difference between the cases χ = χ 0 and χ = χ 3 , except, possibly, for the number of turning points of the secondary loops emanating from the principal curve of the diagram. In particular, all the features described by Theorems 1.1, 1.2 and 1.3 are preserved; among them, the spiraling effect of the principal curve from which the secondary loops emanate as λ decreases.
Naturally, in the case χ = χ 3 , the solution of type 0 on the principal curve cannot exist, as a b is no longer constant in (α, 1−α), however a new clustering phenomenon of points where the solutions change type, represented with small squares on the diagrams, occurs. For example, in the top right plot of Figure 12, the simplest one, following the curve starting at b = 0, we first have solutions of type 1. Then, on the upper half-branch of curve, the type of the solutions becomes five at b = 6.4753, three at b = 6.3486 and, finally, it gets the original value one after b = 5.5997. The plots of Figure 13 (A) represent, from the bottom to the top, solutions on this portion of the diagram with a high number of type-changes, precisely on the upper half-branch for b = 6.4836 (of type 1), for b = 6.4120 (type 5) and for 6.3416 (type 3), respectively.
Another new phenomenon observed in the case χ = χ 3 is that the number of turning points on the secondary loops increases as λ decreases, which is illustrated in the third and fourth plots on the right column of Figure 12. Indeed, first an inflection point on the upper branch of the loop is formed and then two additional turning points arise from it. As a consequence, two other turning points arise on the lower part of the loop, because the asymmetric solutions must occur in pairs of the form (u(x), u(1 − x)). As a consequence, in the last plot of Figure 12, the problem admits at least 8 solutions if 1798.9437 < b < 5303.5426. Five of them, for the special value b = 4300, have been represented in Figure 13 (B). Among them, the two plotted with continuous lines are the ones on the principal curve, and the   Comparing the second row with the first one, one can observe a dramatic change in the bifurcation diagrams, since each of them now consists of two components. Namely, the primary curve, plotted with a continuous (blue) line, and an isola (compact component), plotted with a dashed (red) line. Hence, breaking the symmetry of the weight function provokes the emergence of isolas, instead of loops bifurcating from the primary branch, as it occurs in all the previous (symmetric) situations. The associated breaking of the pitchfork bifurcations has been studied, at a local level, as imperfect bifurcations in, for example, M. Golubitsky and D. Schaeffer [20], R. Gómez-Reñasco and J. López-Gómez [22], and L. Ping, J. Shi and Y. Wang [43].
It should be observed that the isola perturbs from the arc of the principal curve in between the bifurcation points of the loop emerged in case χ = χ 3 and the lower half-arc of this loop if ε < 1 (see the bottom left plot of Figure 14), or the upper half-arc of the loop if ε > 1 (see the bottom right plot of Figure 14). According to these patterns, the unstable manifolds of the solutions along the upper arc of the isola are two-dimensional, while they are one-dimensional on the lower arc, if ε < 1. Similarly, they are two-dimensional on the lower arc and one-dimensional on the upper one if ε > 1. Actually, the theory of Kato [25] can be used to prove that all these dimensions vary analytically with respect to the parameter ε, as the spectra of the linearized equations consist of simple algebraic eigenvalues. It should be also noted that there is a continuous transition among the several points where the solutions change of type along these arcs when ε perturbs from 1.
As a consequence of the continuous dependence, though the topological structure of the bifurcation diagram plotted in the first row of Figure 14 differs from those of the second row, as the number of connected components changes, from a quantitative point of view these diagrams do not differ substantially when ε perturbs from 1. Naturally, these differences are magnified as ε separates away from 1. Indeed, as ε grows the isolas shrink to a single point up to disappear. This striking phenomenon has been depicted in Figure 15, while Table 1 provides us with some relevant quantitative information related to it.  Figure 15. A series of isolas for χ = χ 3,ε with ε = 1.01, 2, 3, 5, 7, 10, 11. As ε grows, they shrink to a single point and disappear.
As far as the primary curves is concerned, the quantitative differences as ε separates away from 1 are not that emphasized, as it becomes apparent by comparing the left plot of Figure 16, which represents the bifurcation diagram for λ = −200 and ε = 11.0084, and the bottom right plot of Figure 14, which gives the corresponding bifurcation diagram for λ = −200 and ε = 1.001. In the former case the isola has almost completely shrunk and the turning point on the principal curve arises at b = 132.1283, while in the latter one the turning point arises at b = 134.0460.
The fact that these isolas shrink to a single point as ε separates away from 1 is far from being in conflict with the multiplicity result of Theorem 1.1, since we can recover the shrunk loops by choosing appropriate (smaller) values of λ. This becomes apparent by comparing the two diagrams of Figure 16. Both of them Table 1. The size of the isolas of Figure 15, for a series of values of ε. were computed for ε = 11.0084, but in the first one we used λ = −200 and in the second one λ = −230. In general, the smaller λ, the larger the isola. Therefore, the expanding tendency of the loops revealed by the numerics of Section 2, in the symmetric case, as λ decreases, is preserved by the isolas arising for general asymmetric χ's. The larger −λ, the wider the isolas. 4. Discretization, numerical codes and computational troubles. To compute the positive solutions of (1) we have coupled a pure spectral method with collocation and a path-following solver. The spectral method uses trigonometric modes. This gives high accuracy at a very low computational cost (see, e.g., Canuto et al. [11]). If N ≥ 1 stands for the number of modes, then, in all our computations we have used N + 2 equidistant collocation points. Set δ = 1/(N + 1), and denote by the collocation points. The number of modes is always chosen in such a way that x i / ∈ {α, 1 − α} for all 0 ≤ i ≤ N , so that, changing the nonlinearity of (1) in the intervals containing α and 1 − α, we can assume, without loss of generality, that the nonlinearity of the differential equation is C 2 . Then, the solutions u(x) of (1) are approximated by the truncated Fourier series and we impose these truncations to satisfy the original boundary value problem at the collocation points, which gives rise to where v = (v 1 , ..., v N ), and we have denoted For general Galerkin approximations, the local convergence as N ↑ ∞ of paths at regular, turning and simple bifurcation points was proven by Brezzi, Rappaz and Raviart [7][8][9]. Incidentally, also convergence for regular singularities of codimension two was proven by López-Gómez, Molina-Meyer and Villareal [40] and by López-Gómez et al. [33]. In all these situations, the local topological structures of the solution curves and/or surfaces for both models, the continuous and the discrete one, are known to be equivalent.
Since u ∈ C 1 [0, 1], its j-th Fourier coefficient v j decays like O(j −1 ) as j ↑ ∞ (see, e.g., Canuto et al. [11]) and, hence, as N ↑ ∞. Due to these features we have used the following criterion to choose the number of modes in all our computations To respect it we needed N = 300 modes in all our numerical computations. The global continuation solvers used to compute the solution curves and the dimensions of the unstable manifolds of all the solutions along them have been adapted from those of Eilbeck [17], Keller [26], López-Gómez [27], Allgower and Georg [2], Crouzeix and Rappaz [16] and López-Gómez et al. [33], where the interested readers are sent for details.
When implementing the available continuation methods in the literature, we found, essentially, two main difficulties. Namely, one must be extremely careful in choosing the shot direction to compute the bifurcated closed loops from the primary curve, and, in addition, one should adopt an appropriate re-scaling procedure to compute automatically all turning points along the primary curve, as they are extremely closed and, hence, the available algorithms in the specialized literature do not work when −λ > 0 is sufficiently large. 4.1. Treatment of the bifurcation directions of the loops. As far as the first difficulty is concerned, it turns out that at the bifurcation points one cannot shoot in an arbitrary non-tangential direction to get the first predicted point on the bifurcated loop. In particular, shooting in the orthogonal direction to the primary curve does not work in our model, because the tangents to the primary curve and to the bifurcated loop are very close to each other in R N +1 at the bifurcation point for large −λ > 0; doing so, the iterates go back to the primary branch again. Consequently, the Method II of Keller [26,Chapter V] does not work in our present situation; this was the method used in most of the global continuations in López-Gómez et al. [33], López-Gómez and Molina-Meyer [35][36][37][38], and Gómez-Reñasco and López-Gómez [22].
Consequently, to compute the bifurcated loops from the primary curve we had to determine very accurately the tangent vectors to the emanated curves at their bifurcation points from the primary curve, as suggested by Keller [26,Method I]. Subsequently, we will detail the way we proceeded. Equation (8) can be expressed as a nonlinear equation Suppose we have already computed the primary curve, say (b(s), v(s)), parameterized by the pseudo-arc length s, around some bifurcation point (b 0 , v 0 ) := (b(0), v(0)). Then, there exists ε > 0 such that and differentiating (9) with respect to s, we are driven to where " · " stands for d/ds. Consequently, Actually, this is the tangent vector to the primary curve at the bifurcation point. It can be easily determined during the computation of the principal curve (b(s), v(s)).
In order to compute the tangent vector to the bifurcated curve, we first need to find a basis of this kernel. Precisely, we want to determine a vector, denoted by Φ, orthogonal to ψ and normalized so that Φ = 1. Thus, setting we should solve BΦ = 0. Unfortunately, in practice, the matrix B is unknown, as an effect of approximation errors. Actually, B is approximated by some matrix A which might be invertible. Therefore, in general AΦ = 0, though, by continuous dependence, AΦ ∼ 0 and A exhibits an eigenvalue perturbed from zero. Therefore, in our numerical experiments, we have taken as Φ the unique φ minimizing A · on the unit surface, i.e., Av .
Note that Φ = φ if A = B. This scheme has shown to be extremely efficient to do all our numerical computations. Up to the best of our knowledge, it goes back to Allgower and Georg [2]. It is well known that φ must be an eigenvector of the symmetric matrix A T A associated to its lowest eigenvalue. Consequently, the inverse power method applied to A T A provides us with an extremely good approximation of φ.
Similarly, one can determineψ ∈ R N such that Using these notations, for any parametrization Y (s), s ∼ 0, of the bifurcated curve with Y (0) = (b 0 , v 0 ), its tangent vector at the bifurcation point is given bẏ  [26,Section 5.26] for the details of the derivation of this formula). As as ε ↓ 0. Naturally, (11) requires much lower computational cost than (10), as it involves only first order derivatives. So, for the numerics of this paper, we have taken the approximatioṅ for an appropriate choice of the auxiliary parameter ε. Our choice was ε = 10 −9 for all the computations.

4.2.
Treatment of the closed turning points. As far as the turning points along the primary branch are concerned, the main difficulty is the exiguous difference between the solutions in each of the half-curves generating the turning point. This causes all the turning points in the bifurcation diagrams to be very closed. Let (b(s), v(s)), s ∈ R, denote the parametrization by pseudo-arc length of the primary curve and pick s 1 < s 2 such that P 1 := (b(s 1 ), v(s 1 )) and P 2 := (b(s 2 ), v(s 2 )) are close to some turning point. Then, the tangent vectors to the primary curve at these points are given by ψ j := (ḃ(s j ),v(s j )), j ∈ {1, 2}, and, in our solvers, the sign of the Euclidean scalar product decides whether P 1 and P 2 are on the same half-curve of the turning point, or not. More precisely, P 1 and P 2 are on the same half-branch if ψ 1 , ψ 2 > 0, while they leave the turning point in between if ψ 1 , ψ 2 < 0.
To compute the tangent vector at s, ψ := (ḃ(s),v(s)), we should solve the equation subject to the normalization condition Solving (13) forv and substituting the result in (14) giveṡ Operatively, one starts by computingḃ from (15), then substitutes the resulting value in (13) and, finally, solves (13) forv. However, we have to choose one of the signs in identity (15). This is the first step in running our path-following solver, depending on whether we want to compute solutions of (1) for increasing values of the bifurcation parameter b, or for decreasing ones. Once fixed the sign ofḃ, it will be kept invariant by the code until a turning point is detected through the negativity of (12). When this occurs, the code changes automatically the sign ofḃ to compute, going backwards, the other half-curve of the turning point. When the solutions are approaching a turning point, the termḃ(s 1 )ḃ(s 2 ) > 0 must converge to zero and, hence, the change of sign of (12) can be detected through the change of sign of the scalar product v(s 1 ),v(s 2 ) . But the converse result is far from true, because v(s 1 ),v(s 2 ) could also change sign at the critical points of the curve, whereḃ(s 1 )ḃ(s 2 ) > 0 is bounded away from zero. Consequently, at the turning points v(s 1 ),v(s 2 ) changes sign andḃ(s 1 )ḃ(s 2 ) ∼ 0. However one can distinguish a turning point from a critical point either by the significant reduction of the magnitude of the continuation step as the solutions are approximating a singular point of F, or through the changes in the dimensions of the unstable manifolds. The most versatile criterion from the practical point of view is the latter, which is the one that we have adopted in this work.
It turns out that, for λ ≤ −800, the turning points are extremely closed. Indeed, the solutions along each of the two half-branches differ from each other by less than 10 −4 . As a consequence, the correction step of the path-following solver might switch from one half-branch to the other one when v(s 1 ),v(s 2 ) < 0 but ψ 1 , ψ 2 > 0. Althoughḃ(s 1 )ḃ(s 2 ) is small, because the continuation step shortened automatically up to order 10 −7 ,ḃ(s 1 )ḃ(s 2 ) is not sufficiently small to provoke the change of sign of (12). In these cases, we know that the solution has switched to the other half-branch by comparing the dimensions of the unstable manifolds of the last two solutions. Such dimensions can change either by one or two unities, according to the relative positions of the last correction with respect to the turning point and the bifurcation points of the secondary loops from the primary curve. Nevertheless, even if at a turning point the quantityḃ(s 1 )ḃ(s 2 ) goes to zero, the fact that it remained large with respect to v(s 1 ),v(s 2 ) when the correction jumped to the other half-branch does not necessarily entail that we are far away from the turning point, because, as a result of the numerics,ḃ(s 1 )ḃ(s 2 ) varies drastically in a neighborhood of the turning point as λ ≤ −800.
Once detected a jump along the primary curve, one must change the sign oḟ b. Otherwise, again by the closeness of the turning point, the algorithm is pushed forward to find new solutions in a small neighborhood of the turning point, where, as a result of the number of solutions available, up to four for sufficiently large −λ, either the corrections might jump to another curve again, or there might be slow convergence of the solver in the correction, causing a severe reduction of the continuation step until reaching the minimal value permitted, which leads to the stop of the algorithm. It should be noted that the turning points and the bifurcation points of the secondary loops from the primary curve do approximate each other as λ ↓ −∞. Consequently, in a neighborhood of a turning point the solver might need to discriminate among four extremely close solutions. As a matter of fact, the path-following solver can provide even with a solution on the bifurcated loop. In these degenerate situations, one should ascertain the symmetry of the solution, besides the dimensions of the unstable manifolds, to determine its precise location. It should be remembered that, according to [41], all solutions along the primary curve are symmetric, while the solutions on the bifurcated loops are asymmetric. Table 2 measures the distance between the value of b where the solution jumps to the other half-curve near the first turning point of the primary branch and the value of b where the solver stops if the sign ofḃ is not changed in (15). We took the latter one as the value of b for the first turning point of the primary curve in all the numerical examples of Section 2. By analyzing Table 2, it is apparent that we are quite close to the turning point when the algorithm jumps to a new branch, confirming thatḃ(s 1 )ḃ(s 2 ) varies dramatically in a neighborhood of a closed turning point.
In many circumstances, when the solution jumps from one branch to another it might cross not only the turning point of the primary branch but also the first bifurcation point of the closed loop emanating from it, which are extremely close for sufficiently large −λ > 0. We know that this occurs because the unstable manifold of the computed solution varied its dimension by 2. In such cases, at a later stage, we refined the computation of the primary curve, going backwards, to approximate as much as possible the bifurcation point of the closed loop. Two possibilities can occur: either we can cross it, or the algorithm jumps to the bifurcated loop. In the second case, we take the first value of b on the loop as the value for the bifurcation point. The larger −λ > 0, the more emphasized the previous trouble, since, according to the numerics, the bifurcation points of the loops approximate the turning points of the primary curve as λ ↓ −∞.
To control automatically the change of sign ofḃ in (15) when a jump along the primary curve occurs, we have used, instead of (12), the weighted scalar product ψ 1 , ψ 2 ε = v(s 1 ),v(s 2 ) + εḃ(s 1 )ḃ(s 2 ) for some appropriate ε > 0, which was chosen according to the size of λ, in order to reduce the weight of the termḃ(s 1 )ḃ(s 2 ). The numerical computations of this paper used the values of ε given in Table 3. For λ < −2000 we were not able to automatize the continuation codes, since we were outside the machine precision range as it is apparent from Table 3.