Multiplicity of nodal solutions in classical non-degenerate logistic equations

: This paper provides a multiplicity result of solutions with one node for a class of (non-degenerate) classical di ff usive logistic equations. Although reminiscent of the multiplicity theorem of L´opez-G´omez and Rabinowitz [1, Cor. 4.1] for the degenerate model, it inherits a completely di ff erent nature; among other conceptual di ff erences, it deals with a di ff erent range of values of the main parameter of the problem. Actually, it is the first existing multiplicity result for nodal solutions of the classical di ff usive logistic equation. To complement our analysis, we have implemented a series of, very illustrative, numerical experiments to show that actually our multiplicity result goes much beyond our analytical predictions. Astonishingly, though the model with a constant weight function can only admit one solution with one interior node, our numerical experiments suggest the existence of non-constant perturbations, arbitrarily close to a constant, with an arbitrarily large number of solutions with one interior node.


Introduction
This paper studies the one-dimensional boundary value problem −u ′′ = λu − a(x)u 3 in (0, 1), u(0) = u(1) = 0, (1.1) where λ ∈ R is a parameter, and a : [0, 1] → [0, ∞), a 0, is a piece-wise continuous function. This problem is said to be non-degenerate if a(x) > 0 for all x ∈ [0, 1], (1. 2) as, in such a case, large positive constants provide us with supersolutions of (1.1) and hence, by the maximum principle, all the solutions (λ, u) with u 0 are bounded. Instead, when [α, β] ⊂ a −1 (0) for some α, β ∈ [0, 1] with α < β, problem (1.1) is said to be degenerate. In these problems, large positive constants fail to be supersolutions, making the analysis of (1.1) substantially harder since large solutions may, and actually do, appear. Degenerate problems have attracted a great deal of attention during the past few decades since the publication of the papers by Brézis and Oswald [2], Ouyang [3] and Fraile et al. [4]. In particular, they have shown to be crucial in developing the theory of metasolutions (see the monograph of López-Gómez [5], as well as the list of references therein). When a −1 (0) is a tiny set with empty interior, (1.1) behaves much like a non-degenerate problem (see, e.g., Daners and López-Gómez [6]).
By adopting the point of view of global bifurcation theory, we see a solution of (1.1) as a pair (λ, u) consisting of a value of λ and a piece-wise C 2 function satisfying the boundary value problem (1.1). Observe that (λ, 0) is a solution for all λ ∈ R; it will be referred to as the trivial solution. The nonlinearity of the differential equation has been chosen to be odd, so that the solutions of (1.1) arise by pairs: (λ, u) and (λ, −u). This considerably simplifies the description of the solution set of (1.1), without reducing the novelty of the findings of this paper.
According to the Cauchy-Lipschitz theorem for second order differential equations, for every solution (λ, u) of (1.1) with u 0 and any z ∈ [0, 1] such that u(z) = 0, necessarily u ′ (z) 0. Therefore, the zeroes of nontrivial solutions, also called nodes in this paper, must be simple. Consequently, u changes sign at any interior node. As a byproduct, the solutions of (1.1) have, at most, finitely many interior nodes, possibly none, as they might be either positive or negative in (0, 1). In this paper, by a nodal solution we mean a solution with n ≥ 1 interior nodes; otherwise, the solution is either positive, or negative. Figure 1 shows some plots of positive, negative and nodal solutions of (1.1) for constant a, whose construction will be described next.  Figure 1. Some nodal solutions of (1.1).
In the simplest case when a(x) is constant, e.g., a ≡ 1, by the symmetries of the differential equation, the one-node solutions satisfying u ′ (0) > 0 can be constructed by simply gluing together the positive solution in (0, 1 2 ) with the negative solution in ( 1 2 , 1). Analogously, the two-node solutions satisfying u ′ (0) > 0 can be constructed by gluing the positive solutions in (0, 1 3 ) and ( 2 3 , 1) with the negative one in ( 1 3 , 2 3 ), and a similar argument can be applied to construct all nodal solutions of (1.1). Indeed, the solutions with n − 1 interior nodes, n ≥ 2, are determined by the positive or negative solutions on appropriate subintervals of the form (a, b), with b − a = 1 n . Such sign-definite solutions exist if, and only if, λ > ( π b−a ) 2 = (nπ) 2 (see, e.g., [7,Theorem 2.1]). When a(x) is constant, this is a very classical result that can be easily established through a shooting argument involving the phase portrait and the symmetries of the differential equation. In the same case, the uniqueness of the nodal solution with u ′ (0) > 0, if it exists, is a direct consequence of the strong maximum principle (see, e.g., Cubillos [8], if necessary). We summarize all these elements in the following well-known result.
Moreover, the set of nontrivial solutions of (1.1) is made up of a sequence of components, C ± n , n ≥ 1, bifurcating from (λ, 0) at λ = λ n , and consisting of solutions with n − 1 interior nodes: C + n denotes the set of solutions (λ, u) with n − 1 nodes such that u ′ (0) > 0, while C − n consists of the solutions with u ′ (0) < 0. By the analyticity of the underlying time maps with respect to λ, along any of these curves, u can be parameterized as an analytic curve with respect to λ, i.e., C ± n are analytic curves for all integer n ≥ 1.
When, instead of being a positive constant, a(x) is assumed to be a general continuous function satisfying (1.2), or a positive piece-wise continuous function, the phase portrait techniques used to derive Theorem 1.1 fail and other techniques, like the topological degree, must be invoked to get the corresponding existence results. In this case, the global results by Rabinowitz [9][10][11] establish that, essentially, the same results are true as soon as (1.2) holds, except for the arc connectedness of the components C ± n . In the setting of [9], C ± n were only shown to be closed and connected subsets, maximal for the inclusion, of the solution set of (1.1), because they were constructed by means of the Leray-Schauder degree. Nevertheless, also in this general case, as a consequence of the strong maximum principle, the positive solution is unique if it exists and, as a byproduct of the implicit function theorem, C + 1 is also a real analytic curve (see, e.g., [4], [5], and [6]). Furthermore, by a recent result of López-Gómez and Sampedro [12], for every n ≥ 2, each of the components C ± n consists of a discrete set of analytic arcs of curve plus a discrete set of branching points, because of the analyticity of the nonlinearity of (1.1). In particular, all these components are arc connected according to the structure theorem of analytic manifolds (see [12] for further details). The reader is sent to Rabinowitz [9][10][11], López-Gómez [13,Chapter 6], and Dancer [14,15] for an overview of the most relevant results in the context of the classical theory.
Although the global bifurcation diagram under condition (1.2) is reminiscent of the one plotted in Figure 2, no uniqueness result of one-node solutions is known for (1.1) under condition (1.2); the uniqueness results of López-Gómez and Rabinowitz [16,17] are only valid in the context of degenerate problems. As a matter of fact, the component C ± 2 might bend backwards providing us with some interval of λ's where (1.1) admit several one-node solutions, or it might even exhibit some secondary, or tertiary, bifurcations, and the existence of some further components of solutions with one interior node cannot be excluded, as those isolas might be provoked by the spatial heterogeneities of a(x). Since the simplest problem of the uniqueness of the one-node solutions under condition (1.2) has remained open since [9][10][11] were published almost fifty years ago, answering this question, either in the affirmative or in the negative, is extremely challenging. This is the main goal of this work. Now, we will introduce the class of weight functions a(x) considered in this paper. First, for any given integer κ ≥ 2, we choose κ subintervals of [0, 1] of the form (α j , β j ), j ∈ {1, ..., κ}, such that • 0 < α 1 < β 1 < α 2 < β 2 < · · · < α κ < β κ < 1, • 0 < h = β j − α j for all j ∈ {1, ..., κ}, i.e., these κ intervals have the same length, h > 0, • α i + β j = 1 if i + j = κ + 1, i.e., (α i , β i ) is the symmetric of (α κ+1−i , β κ+1−i ), with respect to 0.5, for all i ∈ {1, ..., κ}, Then, for each of these choices, we consider the weight functions for every ε ∈ [0, 1], which are symmetric about 1 2 and piece-wise constant. We are considering symmetric weight functions as it is well documented that breaking the symmetry of the problem simply provokes that branches with bifurcation points split into two or more separated components, whose mathematical treatment is more intricate (see, e.g., [18]). In some sense, the symmetries do reorganize the global bifurcation diagrams in a more friendly way. Figure 3 plots two choices of a ε (x) for κ = 4. The first one for ε = 0, the second one for ε = 0.9. Naturally, for the choice ε = 0, problem (1.1) is degenerate with  Figure 3. Two admissible weights a ε (x) with ε = 0 (left) and ε = 0.9 (right).
When κ = 2n for some integer n ≥ 1, we have that β n < 1 2 < α n+1 and, hence, a ε (x) = 1 for all x ∈ (β n , α n+1 ). When, κ = 2n + 1 for some integer n ≥ 1, necessarily, by the symmetry of a ε , one has that α n+1 < 1 2 < β n+1 and hence, a ε (x) = ε for all x ∈ [α n+1 , β n+1 ]. The weight functions a ε (x) defined in (1.3) can be chosen as close as we wish to the constant 1 in [0, 1] by taking ε sufficiently close to 1. Although establishing the precise meaning of 'sufficiently close' might be hard, most of us will agree that ε = 0.9999 is really close to 1. Therefore, for problem (1.1) with a(x) = a 0.9999 (x) one expects, in agreement with Theorem 1.1, a unique one-node solution for some interval of values of λ, regardless the value of κ and α j , β j , j ∈ {1, ..., κ}. However, for the special choice (1.4) with ε = 0.9999, our numerical experiments presented in Section 6 reveal that (1.1) possesses, for some ranges of λ, at least 14 solutions with one interior node: 7 with u ′ (0) > 0, and another 7 with u ′ (0) < 0. Among the first ones, one is symmetric about 1 2 , and the remaining 6 are asymmetric, and 3 among the asymmetric ones are reflection about 1 2 of the remaining 3. Even more astonishingly, our numerical experiments suggest that this multiplicity result holds for really huge intervals of λ's! More precisely, in this paper we are constructing some examples with κ ∈ {2, 3, 4}, for which max such that (1.1) possesses 6, 10 and 14 one-node solutions, respectively, for large intervals of values of λ. Based on our numerical experiments, we conjecture that, whenever a(x) satisfies (1.3), there exists ε 0 > 0 such that, for every ε ∈ (ε 0 , 1), the problem (1.1) possesses (at least) 4κ − 2 one-node solutions for large intervals of λ's. Nevertheless, it is not clear how these branches of solution disappear as ε ↑ 1 to recover the uniqueness that occurs for ε = 1 (see the discussion at the end of Section 4). This is a rather unexpected effect on the global structure of the set of one-node solutions of (1.1), and it is due to the presence of spatial heterogeneities in the weight function a(x). It is striking that one can construct situations with a(x) arbitrarily closed to 1 in the L ∞ -norm for which the problem (1.1) can have an arbitrarily large number of one-node solutions. As a byproduct, the associated global bifurcating diagrams of one-node solutions of (1.1) in the general case when a(x) is not constant have nothing to do with the classical paradigm plotted in Figure 2 (see the global bifurcations diagrams in Sections 4, 5 and 6). More precisely, besides C ± 2 , they might contain up to 2κ − 2 unbounded isolas consisting of asymmetric solutions, and C ± 2 might have one, or more, secondary bifurcations. Throughout this paper, for every p, q ∈ [0, 1] with p < q, and any V ∈ L ∞ (p, q), we will denote by τ n [−D 2 + V; (p, q)], n ≥ 1, the sequence of (real) eigenvalues of the linear eigenvalue problem −φ ′′ + Vφ = τφ in (p, q), φ(p) = φ(q) = 0, ordered so that τ n < τ n+1 for all n ≥ 1. It is well known that each of these eigenvalues is algebraically simple. Moreover, the eigenfunctions associated with τ n [−D 2 + V; (p, q)] have n − 1 interior nodes in (p, q). Thus, for every n ≥ 1, τ n [−D 2 + V; (p, q)] is the unique eigenvalue to which an eigenfunction with n − 1 interior nodes is associated. For notational simplicity, we will denote When a(x) is given by (1.3) with ε = 0, (1.1) becomes a degenerate problem where a −1 (0) has κ ≥ 2 components having the same length h. Assume that h < 1 2 : in such a case, thanks to [7, Theorem 2.1], (1.1) has a positive, or negative, solution, if and only if Moreover, the solution is unique if it exists. Similarly, by [16, Corollary 3.6], (1.1) has a solution with one interior node if and only if λ 2 = (2π) 2 < λ < σ 1 (note that λ 2 < σ 1 , because h < 1 2 ). Furthermore, thanks to [1, Corollary 4.1], there exists η > 0 such that, for every λ ∈ [σ 1 − η, σ 1 ), (1.1) possesses, at least, 4κ − 2 solutions with one interior node.
In this paper, re-elaborating on some of the technical devices of [1], we will show (see Theorem 3.1) that there exists ε 0 > 0 such that, for every ε ∈ (0, ε 0 ], the problem (1.1) with a(x) given by (1.3) possesses at least 2κ − 2 one-node solutions for every λ > σ 1 and at least 4κ − 2 for each λ ∈ (σ 1 , σ 2 ). This result is the first existing multiplicity result of nodal solutions for problem (1.1) in the classical case. Although some ideas are borrowed from [1], our result has a completely different nature, since it deals with the non-degenerate case at a different range of λ's.
Moreover, rather unexpectedly, the numerical results of this paper reveal that, in most of the circumstances, this local multiplicity result is maintained for values of ε very close to 1 and for large intervals of λ. The singular perturbation problem of determining the precise behavior of the global bifurcation diagrams of the one node solutions of (1.1) as ε ↑ 1 remains open in this paper. Our numerical experiments reveal that this problem is extremely intricate.
A crucial aspect of our study is ascertaining the precise evolution of the interior nodes of the onenode solutions of (1.1) as the parameters λ and ε change. Actually, by the uniqueness of the positive solutions of (1.1) on any subinterval of [0, 1], the node, z, of a one-node solution, u, determines it univocally. Thus, z can be identified with u for all practical purposes. By the symmetries of (1.1), when a ≡ 1 it turns out that z = 1 2 for all λ > λ 2 . However, when ε < 1, besides the two symmetric nodal solutions with z = 1 2 , there arise one-node solutions with their respective nodes placed on each of the intervals (α j , β j ), j ∈ {1, ..., κ}.
To conclude this section, we discuss the linearized stability of any one-node solution (λ,ū) of (1.1), which is given by the sign of τ 1 [−D 2 + 3a(x)ū 2 −λ; (0, 1)]. Since and hence, by the monotonicity of the τ n 's with respect to the potential V (see [20], if necessary), we find that Moreover, in the latter case, its unstable manifold is one-dimensional. This paper is organized as follows. Section 2 borrows from López-Gómez and Rabinowitz [1, 16,17] the main concepts and technical devices used in this work. Section 3 establishes the main new analytical result, Theorem 3.1. The rest of the paper discusses a series of numerical experiments to strongly suggest how the local theorems of Sections 2 and 3 are indeed global, in the sense that they still hold true for ε very close to 1 and very large intervals of λ.

Preliminaries. A multiplicity result for ε = 0
Throughout the rest of this paper, we assume that a = a ε has been chosen to satisfy (1.3) with ε ∈ [0, 1]. This problem has been analyzed for the case ε = 0, in López-Gómez and Rabinowitz [1, 16], and López-Gómez, Molina-Meyer and Rabinowitz [7], where, in order to study the solutions of (1.1) with one interior node satisfying u ′ (0) > 0, the following C 1 -map was introduced and w(x; z, λ, ε) is the unique negative solution of provided these solutions exist. Naturally, should they exist, the function x ∈ (z, 1], provides us with a 1-node solution of (1.1) if, and only if, φ(z, λ, ε) = 0. Thus, through this scheme one can identity the 1-node solution u with the position of its node z. In particular, the problem of analyzing the 1-node solutions of (1.1) is equivalent to the problem of analyzing φ −1 (0). When ε > 0, the problems (2.2) and (2.3) are non-degenerate and, hence, v and w are well defined if and only if λ > π z 2 and λ > π 1 − z 2 , respectively. When ε = 0 the situation is more subtle because Thus, the problem is degenerate in the sense that a(x) vanishes somewhere in (0, 1). In such a case, according to [7, Theorem 2.1], whether or not v and w do actually exist depends on the precise location of z. Precisely, Similarly, • when z ≥ β κ , w exists if and only if λ > π 1 − z 2 ; • when z ∈ (α κ , β κ ), w exists if and only if π In particular, this shows that, for every ε ∈ [0, 1], the condition is necessary for the existence of a 1-node solution of (1.1). Note that (z, λ) = 1 2 , λ 2 is the crossing point of the curves λ = ( π z ) 2 and λ = ( π 1−z ) 2 in the (z, λ)-plane. Thus, (2.5) implies that, necessarily, λ > λ 2 . Moreover, in the degenerate case ε = 0, since κ ≥ 2, another necessary condition for the existence of v and w is λ < σ 1 . Thus, when ε = 0, if v and w exist, necessarily (z, λ) ∈ R, where R denotes the set of pairs (z, λ) ∈ (0, 1) × (λ 2 , σ 1 ) satisfying (2.5) that has been depicted in Figure 4. Furthermore, based on the theorem of Crandall and Rabinowitz [19], which explains why (1.1) does have a 1-node solution in this case if, and only if, λ ∈ (λ 2 , σ 1 ), when ε = 0, and if and only if λ > λ 2 , when ε > 0. This methodology goes back to López-Gómez and Rabinowitz [16], where it was developed to characterize the existence of 1-node solutions for the degenerate problem when ε = 0.
where u is a solution of (1.1) with the single interior node z. Since F is continuous, it maps connected subsets in R × C 1 [0, 1] to connected subsets in R 2 . In particular, the component of solutions C + 2 (respectively, C − 2 ) bifurcating from (λ 2 , 0) is transformed into a connected subset denoted by G + 2 (respectively, G − 2 ) bifurcating from ( 1 2 , λ 2 ). Points of the form (z, σ 1 ) may not lie in the range of this map if ε = 0, although they lie in its closure. Figure 4 shows the region R and the component G + 2 in a (real) case where z = 1 2 for all (λ, u) ∈ C + 2 and ε = 0. Figure 4. The region R, for ε = 0, and the sign of φ on its lateral boundaries.
Observe that, although for this particular problem the uniqueness of ℓ 1 can be inferred from the phaseportrait techniques of López-Gómez and Rabinowitz [1], it is a direct consequence from the multidimensional uniqueness theorems of López-Gómez and Maire [21] and López-Gómez, Maire and Véron [22].
The limiting behavior of these two functions depends on the precise location of z in (0, 1). The next preliminary result shows that, whenever z ≥ β j , j ≥ 1, the function v(·; z, λ, ε) blows up in (α i , β i ) for all i ∈ {1, ..., j} as ε ↓ 0. It extends some previous findings of [24] in a different context.
Naturally, the negative solution w ε (x) := w(x; z, λ, ε) satisfies similar properties as v ε when ε ↓ 0. Thus, to obtain the proof of Theorem 3.1 from Proposition 3.4 and Remark 3.5, one can argue as follows.
In agreement with Theorem 2.2, (1.1) has three 1-node solutions with u ′ (0) > 0 for each λ < σ 1 sufficiently close to σ 1 . One of them satisfies z = 0.5 and is odd about 0.5. Numerically, we see that the other 1-node solutions satisfy z ∼ 0. 35  We now describe the stability of the solutions, which determines their local attractive or repelling character. As noticed in Section 1, we already know that the linear stability of any 1-node solution is determined by the sign of the lowest eigenvalue of its linearized equation and that, if a 1-node solution is unstable, then its unstable manifold is one-dimensional. According to [19], the branch of nodal solutions with z = 0.5 bifurcates from (λ, u) = (λ, 0) at λ = λ 2 . Thanks to the exchange stability principle [29], it consists of unstable solutions until λ reaches λ b . At such a point, the solutions with z = 0.5 become stable for any further value of λ. Once again, by the exchange stability principle, 3), respectively, as λ ↑ σ 1 . In Figure 5, as in the remaining bifurcation diagrams of this paper, the dashed curves represent unstable solutions, while the continuous ones consist of stable solutions. Figure 6 shows a series of solutions along each of these curves for a series of increasing values of λ. In the first picture we have plotted the nodal solutions with z = 0.5. It is apparent that, as λ increases, they become much larger in  In the right picture of Figure 6 we have plotted a series of solutions on the right bifurcated branch of Figure 5. By the internal symmetries of (1.1) for the choice (1.3), with (4.1), they are reflections about 0.5 of the ones plotted in the central picture. Thus, their nodes stabilize to 0.65, and they blow-up to +∞ in [0.3, 0.4] as λ ↑ σ 1 .
As expected from Theorem 3.1(a), for sufficiently small ε > 0, some of the previous nodal solutions are defined for all values of λ > σ 1 . Surprisingly, our numerical computations show that this is true also for ε close to 1. For example, Figure 7 shows the corresponding global bifurcation diagram for ε = 0.9999. Essentially, the solutions behave in the same way as in the special case when ε = 0, except for the relevant fact that they are defined for all values of λ > σ 1 for which we performed the computations. Figure 8 shows a series of plots of the nodal solutions along each of the three branches of Figure 7. As before, the left picture superimposes the plots of a series of solutions along the central branch, with z = 0.5, the central picture plots a series of solutions on the left secondary branch, whose nodes approach 0.35 as λ increases, and the third one superimposes the plots of a series of solutions on the right secondary branch, whose nodes approach 0.65 as λ grows. The reader should compare the huge qualitative differences between the growth of the solutions plotted in Figures 6 and 8. Although the nodes behave in a rather similar way, approaching either 0.5, or the midpoints of [0.3, 0.4] or [0.6, 0.7], the solutions are much flatter in the present situation, where a 0.9999 is very close to the constant 1. Therefore, although Theorem 3.1 only guarantees, for sufficiently small ε > 0, the existence of one solution for every λ > σ 1 and of three if λ ∈ (σ 1 , σ 2 ), this local result might actually have a global character in the sense that, according to our numerical experiments, there are very large ranges of λ and ε for which (1.1) still has three nodal solutions with one interior node.
And, what is even more striking, is that, according to our numerical simulations, the same structure is maintained at least up to ε = 0.99999999, though in such a case the bifurcation point from the central branch grew up to λ b = 435.0634. This causes a great perplexity to us, as problem (1.1) with a ≡ 1 has a unique solution with one interior node at z = 0.5, whereas the nodes of the solutions along the secondary branches for ε = 0.99999999 still approach 0.35 and 0.65 as λ increases for really huge intervals of λ. The extremely challenging problem of ascertaining how these nodal solutions disappear as ε ↑ 1, and only one solution remains for ε = 1, is left open in this paper.
We can think of two possibilities. The first one is that the secondary branches actually form a loop bifurcating from the central branch at an another (huge) value of λ, and that such a loop shrinks to a single point up to disappear as ε ↑ 1, analogously to what happens in Figure 14. This global structure of the bifurcation diagram could be established, for example, if one is able to adapt to this situation the arguments provided by Dancer and López-Gómez in [30, Theorem 3.1]. The second possibility is that the bifurcation point λ b (ε) converges to ∞ as ε ↑ 1. Actually, analyzing this singular perturbation problem is a serious challenge which needs to be tackled if one wishes to understand how local heterogeneities can have global effects even at large parameter scales. The challenge arises also from the numerical point of view, as, for sufficiently large λ, we are working close to the computational limits of the computer.
In this section we present the results of our numerical computation related to a weight a ε as in (1.3) with κ = 3 and In particular, (4.2) still holds true. Figure 9 shows the global bifurcation diagram of 1-node solutions that we have computed for ε = 0. In agreement with Theorem 2.2, problem (1.1) has 5 nodal solutions (λ, u) with u ′ (0) > 0 for λ < σ 1 sufficiently close to σ 1 . Nevertheless, our computations show that there is some intermediate range of values of λ where the problem has up to seven solutions with u ′ (0) > 0. By simply having a glance at Figure 9, it is easily realized that the set of nodal solutions of (1.1) with u ′ (0) > 0 consists of three components. Namely, the central one G + 2 , plus two additional components: one on the left of G + 2 , say H + 2,ℓ , and another one on its right, H + 2,r , which are isolated with respect to the trivial solution (λ, u) = (λ, 0). Thus, for every λ ∈ (λ t , λ b,2 ) problem (1.1) has, at least, 14 solutions, (λ, u): 7 among them satisfy u ′ (0) > 0, and their opposite, with u ′ (0) < 0. The solutions with z = 0.5 are unstable for λ ∈ (λ 2 , λ b,1 ), stable for λ ∈ (λ b,1 , λ b,2 ) and unstable for λ ∈ (λ 2,b , σ 1 ), in agreement with the exchange stability principle, [29]. Figure 10 shows the plots of a series of 1-node solutions with z = 0.5 along G + 2 . The solutions on the first picture were computed for λ close to λ 2 = (2π) 2 , while the solutions on the second picture where computed for a range of λ's closer to σ 1 . These solutions blow up to +∞ in [0.2, 0.3], and to −∞ in [0.7, 0.8], as λ ↑ σ 1 , whereas they stabilize on the complement of these two intervals to the appropriate large solutions, either positive, negative, or nodal.   Finally, Figure 13 shows the plots of a series of nodal solutions on the component H + 2,r for the same range of values of λ used in Figure 12. The profiles are similar to those computed along H + 2,ℓ , except for the fact that now the nodes approach the values 0.625 and 0.75, respectively, as λ ↑ σ 1 . In particular, these solutions blow up to +∞ in [0.2, 0.3] ∪ [0.45, 0.55] as λ ↑ σ 1 . Moreover, the stable solutions, those with z ∼ 0.625, also blow up to −∞ in [0, 7, 0.8] as λ ↑ σ 1 .
This precise point-wise behavior of the nodal solutions as λ ↑ σ 1 agrees with the results of Theorems 2.1 and 2.2. Actually, such an information has been crucial to construct the global bifurcation diagram plotted in Figure 9, as it will become apparent in a few moments. To compute the component G + (or, equivalently, C + 2 ), we numerically solved the bifurcation problem from u = 0 at λ = λ 2 and then continued the bifurcated curve, detecting and marking the two bifurcation points along it. We stopped when we reached a value of λ sufficiently close to σ 1 . Once computed the whole branch in this way, we went back to the secondary bifurcation points to complete the computation of the closed loop.
In the global bifurcation diagram for ε = 0.1, the loop of solutions on G + 2 bifurcates from the solution with z = 0.5 at λ b,1 = 69.5268 and λ b,2 = 193.3319, while the turning points of the components H + 2,ℓ and H + 2,r occur at λ t = 173.5823. In the global bifurcation diagram for ε = 0.7, these values changed to λ b,1 = 110.7253, λ b,2 = 160.4316, and λ t = 214.8051. According to our numerical experiments, the closed loop bifurcating from the nodal solutions with z = 0.5, is persistent for all values of ε until some critical value ε c ∈ (0.77, 0.78) where it shrinks to a single point before disappearing for all further values of ε for which we have computed the global bifurcation diagram. Once the closed loop disappears, all the solutions of the central branch are unstable. In the case ε = 0.8, λ t = 230.2345. Finally, λ t = 1.1068 × 10 3 when ε = 0.9999. Thus, it is apparent that, although λ t increases as the secondary bifurcation parameter ε approaches 1, the increasing rate is not as high as one might expect. Anyway, observe that the structure of the bifurcation diagrams is maintained for all the larger values of ε for which we have computed them, going much beyond the local perturbation result given by Theorem 3.1.
6. An example with κ = 4 In this section we will discuss the model (1.1) for the choice (1.3) with κ = 4, ε ∈ [0, 1) and (1.4). According to Theorems 2.2 and 3.1 the associated problems admit, at least, 14 solutions with one interior node, for λ in a left neighborhood of σ 1 when ε = 0 and for λ ∈ (σ 1 , σ 2 ) and ε > 0 sufficiently small, respectively. Moreover, for λ ≥ σ 2 and ε > 0 sufficiently small, Theorem 3.1 guarantees the existence of, at least, 6 solutions with one interior node. In Figure 15 we have plotted the global bifurcation diagram that we have computed for ε = 0. In agreement with the result of Theorem 2.2, for the present choice of a(x), problem (1.1) possesses, at least, seven solutions (λ, u) with one interior node and u ′ (0) > 0. By the symmetry of the problem, the pairs (λ, −u) provide us with other seven (different) solutions.
As in the case analyzed in Section 5, the solution set consists of three components. namely, the component, G + 2 (or, equivalently, C + 2 ), bifurcating from u = 0 at λ = λ 2 , plus two additional supercritical folds, H + 2,ℓ and H + 2,r . The component G + 2 consists of the "symmetric"solutions with z = 0.5 plus a secondary curve of solutions bifurcating supercritically from the primary curve at λ b = 74.8175 and filled in by linearly unstable solutions, in full agreement with the exchange stability principle, [29]. As the turning points of the components H + 2,ℓ and H + 2,r occur at λ t = 433.9062 > λ b , problem (1.1) has -limited to solutions with one interior node and u ′ (0) > 0 -, at least, one solution (λ, u) for each λ ∈ (λ 2 , λ b ], three solutions for each λ ∈ (λ b , λ t ), five at λ = λ t and seven for each λ ∈ (λ t , σ 1 ). The nodes of these solutions, as well as their point-wise behaviors as λ ↑ σ 1 , obey the general patterns described in Section 2. Therefore, the solution blows up to ±∞ as λ ↑ σ 1 in every component of a −1 0 (0), except the one containing the node z. Precisely, if u ′ (0) > 0 and the component of a −1 0 (0) completely lies on the left of z, the solution blows up to +∞, while it approximates −∞ if the component completely lies on the right of z.
As suggested by Theorem 3.1, the structure of the global bifurcation diagram plotted in Figure 16 is maintained for ε > 0, ε ∼ 0, at least for λ ∈ (σ 1 , σ 2 ). Nevertheless, our numerical experiments show that the same global structure persists for all further values of ε ∈ (0, 1) for which we have computed it, and the full multiplicity result given by Theorem 3.1(b) also holds for λ ≥ σ 2 . Figure 16 shows a series of global bifurcation diagrams computed for ε = 0.1, ε = 0.5, ε = 0.9 and ε = 0.9999, confirming that our multiplicity results, of a local nature, are actually valid, at least, for all ε ∈ (0, 0.9999].  Table 1 provides us with the values of ε, λ b and λ t in each of the global bifurcation diagrams plotted in Figure 16.
As in all cases analyzed in the previous sections, the analysis of the fine behavior of the global bifurcation diagram as ε ↑ 1 remains an open problem in this paper.
To conclude, we remark that, in general, the number of components of the solution set of (1.1) can Table 1. The values of ε, λ b (ε) and λ t (ε) in Figure 16. vary with the distribution of the intervals where a ε (x) = ε, as well as with the value of the parameter ε.
In particular, some of our numerical simulations (not presented here) suggest that there are examples where, besides the components G + 2 , H + 2,ℓ and H + 2,r , the model exhibits some additional "isola" also containing nodal solutions. A complete description of these more sophisticated aspects will appear elsewhere.
This method gives high accuracy at a rather reasonable computational cost (see, e.g., Canuto, Hussaini, Quarteroni and Zang [36]). In particular, the pseudo-spectral method is very efficient and versatile for choosing the shooting direction from the trivial solution in order to compute small classical solutions, as it provides us with the true values for the first N bifurcation points from the trivial solution, where N is the number of internal mesh points (see Eilbeck [37]).
However, in this paper, we have preferred to use a centered finite difference scheme with equidistant points in order to discretize (1.1), because it is more adequate for solutions with large gradients, like ours (cf. the right plot of Figure 10). This method runs much faster when computing global solution branches in the bifurcation diagrams than the pseudo-spectral method but gives only rough approximations of the true bifurcation values from the trivial solution. To compensate for this inconvenient, the finite difference method requires a substantially higher number of internal mesh points N in comparison to the pseudo-spectral method. Indeed, the k-th eigenvalue (k ≤ N) of the tridiagonal Toeplitz matrix corresponding to the discretization of the operator − d 2 dx 2 is given bỹ λ k =λ k (N) = 2(N + 1) 2 1 + cos N + 1 − k N + 1 π (see, e.g., [38, Section 2]), and it converges to the exact value λ k = (kπ) 2 as N → +∞. For our computations, we chose N = 501, so that λ 2 −λ 2 < 5.2 · 10 −4 .
For general Galerkin approximations, the local convergence of the solution paths at regular, turning and simple bifurcation points was proven by Brezzi, Rappaz and Raviart in [39][40][41], and by López-Gómez, Molina-Meyer and Villareal [42] and López-Gómez, Eilbeck, Duncan and Molina-Meyer in [32] for codimension two singularities in the context of systems. Such results rely on the approximate implicit function theorem (see, e.g., López-Gómez, [43,Theorem 3.2]). In all these situations, the local structure of the solution sets for the continuous and the discrete models are equivalent provided that a sufficiently fine discretization is performed. The global continuation solvers used to compute the solution curves of this paper, as well as the dimensions of the unstable manifolds of the solutions filling them, have been built from the theory of Allgower and Georg [44], Crouzeix and Rappaz [45], Eilbeck [37], Keller [46], López-Gómez [43] and López-Gómez, Eilbeck, Duncan and Molina-Meyer [32].
In order to compute the component C + 2 (or, equivalently, G + 2 ), we started from a value of λ bigger than, but close to the bifurcation point λ 2 and used the initial iterate u 0 = sin(2πx), since, there, it provides us with a first order approximation of the shape of the solutions bifurcating from (λ 2 , 0), and our correction algorithm based on Newton's method easily converged to a solution of problem (1.1). Then, by performing a global continuation with the Keller-Yang algorithm [47], and detecting whether secondary bifurcation points were present, we were able to obtain the global structure of C + 2 . As for the other components, namely the global supercritical folds, Theorem 2.1 allows us to know the shape of the solutions, in particular the position of their nodes for λ < σ 1 , λ ∼ σ 1 , and we used this information to construct adequate predictions to compute globally these components. For instance, when we computed the component containing the nodal solutions with one interior node z close to 0.25 in Figure 9, we chose as the initial iterate, u 0 , for the underlying Newton method used to compute a first nodal solution on the curve, some bounded, not necessarily continuous function u 0 which is positive in [0, 0.25) and negative in (0. 25,1] or, at least, whose node is close enough to 0.25. In addition, this initial iterate u 0 must be large enough, otherwise the Newton method converges to the trivial solution (λ, 0). An a posteriori analysis seems to indicate that u 0 is adequate provided its discrete L 2 norm is greater than the discrete L 2 norm of the solution we are trying to calculate. Once located this first point on the component, our path-following code provided us with the entire component using the Keller and Yang algorithm [47] to transform the supercritical turning points of these folds, which are singular points for the original system, into regular points for an augmented system, as explained by López-Gómez in [43, p. 206].
The complexity of the computed bifurcation diagrams, due to the existence of interior and boundary layers for the nodal solutions, required an extremely careful control of all the continuation steps in our codes for values of λ close to σ 1 when ε ∼ 0, as well as for all λ > λ 2 = (2π) 2 when ε ∼ 1. This explains why the existing commercial packages, such as AUTO-07P, are of no utility when dealing with differential equations in the presence of spatial heterogeneities. Actually, Doedel and Oldeman already emphasized on page 185 of the AUTO-07P manual available at https://depts.washington.edu/bdecon/workshop2012/auto-tutorial/documentation that "given the non-adaptive spatial discretization, the computational procedure here is not appropriate for PDEs with solutions that rapidly vary in space, and care must be taken to recognize spurious solutions and bifurcations." Indeed, this is one of the main problems that we found in our numerical experiments, because, for