Multiplicity of nodal solutions in classical non-degenerate logistic equations

Dedicated to E. N. Dancer on the occasion of his 75th birthday, with profound respect for his mathematical contributions and our very best wishes for his personal files.

1.   Introduction

This paper studies the one-dimensional boundary value problem

where ${\lambda}\in{{\mathbb R}}$ is a parameter, and $a: [0, 1]\to [0, \infty)$, $a\neq 0$, is a piece-wise continuous function. This problem is said to be non-degenerate if

as, in such a case, large positive constants provide us with supersolutions of (1.1) and hence, by the maximum principle, all the solutions $({\lambda}, u)$ with $u\neq 0$ are bounded. Instead, when $[{\alpha}, {\beta}]\subset a^{-1}(0)$ for some ${\alpha}, {\beta}\in [0, 1]$ with ${\alpha} < {\beta}$, 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 $({\lambda}, u)$ consisting of a value of ${\lambda}$ and a piece-wise $\mathcal {C}^2$ function satisfying the boundary value problem (1.1). Observe that $({\lambda}, 0)$ is a solution for all ${\lambda}\in{{\mathbb 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: $({\lambda}, u)$ and $({\lambda}, -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 $({\lambda}, u)$ of (1.1) with $u\neq 0$ and any $z\in [0, 1]$ such that $u(z) = 0$, necessarily $u'(z)\neq 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\geq 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.

In the simplest case when $a(x)$ is constant, e.g., $a\equiv 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, \frac{1}{2})$ with the negative solution in $(\frac{1}{2}, 1)$. Analogously, the two-node solutions satisfying $u'(0) > 0$ can be constructed by gluing the positive solutions in $(0, \frac{1}{3})$ and $(\frac{2}{3}, 1)$ with the negative one in $(\frac{1}{3}, \frac{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\geq 2$, are determined by the positive or negative solutions on appropriate subintervals of the form $(a, b)$, with $b-a = \frac{1}{n}$. Such sign-definite solutions exist if, and only if, ${\lambda} > (\frac{\pi}{b-a})^2 = (n\pi)^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.

Theorem 1.1. Suppose $a\equiv 1$. Then,

${\rm (a)}$ problem $(1.1)$ has a positive solution if, and only if, ${\lambda} > {\lambda}_1 = : \pi^2$. Moreover, it is symmetric about $x = \frac{1}{2}$;

${\rm (b)}$ for every integer $n\geq 2$, problem $(1.1)$ has a solution with $n-1$ (interior) nodes in $(0, 1)$, if, and only if, ${\lambda} > {\lambda}_n = : (n\pi)^2$. Moreover, in such case, its nodes must be equidistant;

${\rm (c)}$ for every integer $n\geq 1$ and ${\lambda} > (n\pi)^2$, there exists a unique solution of $(1.1)$ with $n-1$ interior nodes such that $u'(0) > 0$.

Figure 2 shows the global bifurcation diagram of (1.1) in the special case when $a\equiv 1$. For any solution $({\lambda}, u)$, we are plotting the value of the parameter ${\lambda}$, in abscissas, versus the value $\max u$ if $u'(0) > 0$, or $\min u$ if $u'(0) < 0$, in ordinates. Thus, each point on every plotted curve represents a different solution of (1.1).

According to Figure 2, $({\lambda}, 0)$ is the unique solution of (1.1) if ${\lambda}\leq {\lambda}_1$, (1.1) has a positive solution and a negative solution for each ${\lambda}\in ({\lambda}_1, {\lambda}_2]$, a positive solution, a negative solution, a one-node solution with $u'(0) > 0$ and another one with $u'(0) < 0$ if ${\lambda} \in ({\lambda}_2, {\lambda}_3]$, and so on… Thus, whenever ${\lambda}$ increases and crosses a value of ${\lambda}_n$, the number of solutions of (1.1) increases by two.

Moreover, the set of nontrivial solutions of (1.1) is made up of a sequence of components, $\mathcal{C}_n^\pm$, $n\geq 1$, bifurcating from $({\lambda}, 0)$ at ${\lambda} = {\lambda}_n$, and consisting of solutions with $n-1$ interior nodes: $\mathcal {C}^+_n$ denotes the set of solutions $({\lambda}, u)$ with $n-1$ nodes such that $u'(0) > 0$, while $\mathcal {C}^-_n$ consists of the solutions with $u'(0) < 0$. By the analyticity of the underlying time maps with respect to ${\lambda}$, along any of these curves, $u$ can be parameterized as an analytic curve with respect to ${\lambda}$, i.e., $\mathcal {C}_n^\pm$ are analytic curves for all integer $n\geq 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 $\mathcal {C}_n^\pm$. In the setting of ^[9], $\mathcal {C}_n^\pm$ 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, $\mathcal {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\geq 2$, each of the components $\mathcal {C}_n^\pm$ 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 $\mathcal {C}_2^\pm$ might bend backwards providing us with some interval of ${\lambda}$'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 $\kappa\geq 2$, we choose $\kappa$ subintervals of $[0, 1]$ of the form $({\alpha}_j, {\beta}_j)$, $j\in\{1, ..., \kappa\}$, such that

● $0 < {\alpha}_1 < {\beta}_1 < {\alpha}_2 < {\beta}_2 < \cdots < {\alpha}_\kappa < {\beta}_\kappa < 1$,

● $0 < h = {\beta}_j-{\alpha}_j$ for all $j\in\{1, ..., \kappa\}$, i.e., these $\kappa$ intervals have the same length, $h > 0$,

● ${\alpha}_i+{\beta}_j = 1$ if $i+j = \kappa+1$, i.e., $({\alpha}_i, {\beta}_i)$ is the symmetric of $({\alpha}_{\kappa+1-i}, {\beta}_{\kappa+1-i})$, with respect to $0.5$, for all $i\in \{1, ..., \kappa\}$,

● $h < {\alpha}_1 = 1-{\beta}_\kappa$.

Then, for each of these choices, we consider the weight functions

for every ${\varepsilon} \in [0, 1]$, which are symmetric about $\tfrac{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_{{\varepsilon}}(x)$ for $\kappa = 4$. The first one for ${\varepsilon} = 0$, the second one for ${\varepsilon} = 0.9$. Naturally, for the choice ${\varepsilon} = 0$, problem (1.1) is degenerate with

Thus, $h = 0.1$ and

When $\kappa = 2n$ for some integer $n\geq 1$, we have that ${\beta}_n < \tfrac{1}{2} < {\alpha}_{n+1}$ and, hence, $a_{\varepsilon}(x) = 1$ for all $x\in ({\beta}_n, {\alpha}_{n+1})$. When, $\kappa = 2n+1$ for some integer $n\geq 1$, necessarily, by the symmetry of $a_{\varepsilon}$, one has that ${\alpha}_{n+1} < \tfrac{1}{2} < {\beta}_{n+1}$ and hence, $a_{\varepsilon}(x) = {\varepsilon}$ for all $x\in [{\alpha}_{n+1}, {\beta}_{n+1}]$.

The weight functions $a_{\varepsilon}(x)$ defined in (1.3) can be chosen as close as we wish to the constant $1$ in $[0, 1]$ by taking ${\varepsilon}$ sufficiently close to $1$. Although establishing the precise meaning of "sufficiently close" might be hard, most of us will agree that ${\varepsilon} = 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 ${\lambda}$, regardless the value of $\kappa$ and ${\alpha}_j$, ${\beta}_j$, $j\in\{1, ..., \kappa\}$. However, for the special choice (1.4) with ${\varepsilon} = 0.9999$, our numerical experiments presented in Section 6 reveal that (1.1) possesses, for some ranges of $\lambda$, 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 $\tfrac{1}{2}$, and the remaining $6$ are asymmetric, and $3$ among the asymmetric ones are reflection about $\tfrac{1}{2}$ of the remaining $3$. Even more astonishingly, our numerical experiments suggest that this multiplicity result holds for really huge intervals of ${\lambda}$'s! More precisely, in this paper we are constructing some examples with $\kappa\in\{2, 3, 4\}$, for which

such that (1.1) possesses $6$, $10$ and $14$ one-node solutions, respectively, for large intervals of values of ${\lambda}$. Based on our numerical experiments, we conjecture that, whenever $a(x)$ satisfies (1.3), there exists ${\varepsilon}_0 > 0$ such that, for every ${\varepsilon}\in ({\varepsilon}_0, 1)$, the problem (1.1) possesses (at least) $4\kappa-2$ one-node solutions for large intervals of ${\lambda}$'s. Nevertheless, it is not clear how these branches of solution disappear as ${\varepsilon}{\uparrow} 1$ to recover the uniqueness that occurs for ${\varepsilon} = 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^\infty$-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 $\mathcal {C}_2^\pm$, they might contain up to $2\kappa-2$ unbounded isolas consisting of asymmetric solutions, and $\mathcal {C}_2^\pm$ might have one, or more, secondary bifurcations.

Throughout this paper, for every $p, q\in [0, 1]$ with $p < q$, and any $V\in L^\infty(p, q)$, we will denote by $\tau_n[-D^2+V; (p, q)]$, $n\geq 1$, the sequence of (real) eigenvalues of the linear eigenvalue problem

ordered so that $\tau_n < \tau_{n+1}$ for all $n\geq 1$. It is well known that each of these eigenvalues is algebraically simple. Moreover, the eigenfunctions associated with $\tau_n[-D^2+V; (p, q)]$ have $n-1$ interior nodes in $(p, q)$. Thus, for every $n\geq 1$, $\tau_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 ${\varepsilon} = 0$, (1.1) becomes a degenerate problem where $a^{-1}(0)$ has $\kappa\geq 2$ components having the same length $h$. Assume that $h < \tfrac{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

(note that ${\lambda}_2 < {\sigma}_1$, because $h < \tfrac{1}{2}$). Furthermore, thanks to [1,Corollary 4.1], there exists $\eta > 0$ such that, for every ${\lambda}\in [{\sigma}_1-\eta, {\sigma}_1)$, (1.1) possesses, at least, $4\kappa-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 ${\varepsilon}_0 > 0$ such that, for every ${\varepsilon}\in (0, {\varepsilon}_0]$, the problem (1.1) with $a(x)$ given by (1.3) possesses at least $2\kappa-2$ one-node solutions for every ${\lambda} > \sigma_1$ and at least $4\kappa-2$ for each ${\lambda}\in(\sigma_1, \sigma_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 ${\lambda}$'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 ${\varepsilon}$ very close to 1 and for large intervals of ${\lambda}$. The singular perturbation problem of determining the precise behavior of the global bifurcation diagrams of the one node solutions of (1.1) as ${\varepsilon}{\uparrow} 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 one-node solutions of (1.1) as the parameters ${\lambda}$ and ${\varepsilon}$ 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\equiv 1$ it turns out that $z = \tfrac{1}{2}$ for all ${\lambda} > {\lambda}_2$. However, when ${\varepsilon} < 1$, besides the two symmetric nodal solutions with $z = \tfrac{1}{2}$, there arise one-node solutions with their respective nodes placed on each of the intervals $({\alpha}_j, {\beta}_j)$, $j\in\{1, ..., \kappa\}$.

To conclude this section, we discuss the linearized stability of any one-node solution $(\bar{\lambda}, \bar u)$ of (1.1), which is given by the sign of $\tau_1[-D^2+3a(x)\bar u^2-\bar {\lambda}; (0, 1)]$. Since

it becomes apparent that

and hence, by the monotonicity of the $\tau_n$'s with respect to the potential $V$ (see ^[20], if necessary), we find that

Therefore, $(\bar{\lambda}, \bar u)$ is linearly stable if

while it is linearly unstable if

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 ${\varepsilon}$ very close to $1$ and very large intervals of ${\lambda}$.

2.   Preliminaries. A multiplicity result for ${\varepsilon}=0$

Throughout the rest of this paper, we assume that $a = a_{\varepsilon}$ has been chosen to satisfy (1.3) with ${\varepsilon}\in [0, 1]$. This problem has been analyzed for the case ${\varepsilon} = 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 $\mathcal {C}^1$-map was introduced

where $v(x; z, {\lambda}, {\varepsilon})$ is the unique positive solution of

and $w(x; z, {\lambda}, {\varepsilon})$ is the unique negative solution of

provided these solutions exist. Naturally, should they exist, the function

provides us with a 1-node solution of (1.1) if, and only if, ${\varphi}(z, {\lambda}, {\varepsilon}) = 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 ${\varphi}^{-1}(0)$.

When ${\varepsilon} > 0$, the problems (2.2) and (2.3) are non-degenerate and, hence, $v$ and $w$ are well defined if and only if

respectively. When ${\varepsilon} = 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,

In particular, this shows that, for every ${\varepsilon}\in[0, 1]$, the condition

is necessary for the existence of a 1-node solution of (1.1). Note that $(z, {\lambda}) = \left(\tfrac{1}{2}, {\lambda}_2\right)$ is the crossing point of the curves ${\lambda} = (\tfrac{\pi}{z})^2$ and ${\lambda} = (\tfrac{\pi}{1-z})^2$ in the $(z, {\lambda})$-plane. Thus, (2.5) implies that, necessarily, ${\lambda} > {\lambda}_2$. Moreover, in the degenerate case ${\varepsilon} = 0$, since $\kappa\geq 2$, another necessary condition for the existence of $v$ and $w$ is ${\lambda} < {\sigma}_1$. Thus, when ${\varepsilon} = 0$, if $v$ and $w$ exist, necessarily $(z, {\lambda})\in \mathcal {R}$, where $\mathcal {R}$ denotes the set of pairs $(z, {\lambda})\in (0, 1)\times ({\lambda}_2, {\sigma}_1)$ satisfying (2.5) that has been depicted in Figure 4. Furthermore, based on the theorem of Crandall and Rabinowitz ^[19],

Therefore, for every ${\varepsilon}\in [0, 1]$,

which explains why (1.1) does have a 1-node solution in this case if, and only if, ${\lambda}\in ({\lambda}_2, {\sigma}_1)$, when ${\varepsilon} = 0$, and if and only if ${\lambda} > {\lambda}_2$, when ${\varepsilon} > 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 ${\varepsilon} = 0$.

For the special choice (1.3) with ${\varepsilon}\in[0, 1]$, as $a(x)$ is symmetric about $x = \tfrac{1}{2}$, it is apparent that

and hence ${\varphi}(\tfrac{1}{2}, {\lambda}, {\varepsilon}) = 0$ for all ${\lambda} \in ({\lambda}_2, {\sigma}_1)$ if ${\varepsilon} = 0$, and for all ${\lambda} > {\lambda}_2$ if ${\varepsilon} > 0$. In particular, (1.1) has a unique 1-node solution with $u'(0) > 0$ and $z = \tfrac{1}{2}$.

Let $\mathcal {F}$ denote the map $({\lambda}, u)\mapsto (z, {\lambda})$, where $u$ is a solution of (1.1) with the single interior node $z$. Since $\mathcal {F}$ is continuous, it maps connected subsets in ${{\mathbb R}}\times \mathcal {C}^1[0, 1]$ to connected subsets in ${{\mathbb R}}^2$. In particular, the component of solutions $\mathcal {C}_2^+$ (respectively, $\mathcal {C}_2^-$) bifurcating from $({\lambda}_2, 0)$ is transformed into a connected subset denoted by $\mathcal {G}_2^+$ (respectively, $\mathcal {G}_2^-$) bifurcating from $(\tfrac{1}{2}, {\lambda}_2)$. Points of the form $(z, {\sigma}_1)$ may not lie in the range of this map if ${\varepsilon} = 0$, although they lie in its closure. Figure 4 shows the region $\mathcal {R}$ and the component $\mathcal {G}_2^+$ in a (real) case where $z = \tfrac{1}{2}$ for all $({\lambda}, u)\in \mathcal {C}_2^+$ and ${\varepsilon} = 0$.

Although, according to (2.4), on the top boundary of $\mathcal {R}$ at level ${\lambda} = {\sigma}_1$ the solution $v(x; z, {\sigma}_1, 0)$ is not defined if $z\geq {\beta}_1$, we can extend the definition of $v(x; z, {\sigma}_1, 0)$ to cover the case $z\geq {\beta}_1$ by setting

Indeed, since $v(x; z, {\lambda}, 0)$ is increasing with respect to ${\lambda}$ for every $x\in (0, z)$, the limit (2.6) is well-defined. Moreover, with the ingredients of the proofs of [7,Theorem 2.3], [1,Theorem 2.2] and [1,Theorem 2.6], the next properties hold.

Suppose $z\in (0, {\beta}_1)$. Then, $v(x; z, {\sigma}_1, 0)$ is the unique positive solution of (2.2) at ${\lambda} = {\sigma}_1$ and ${\varepsilon} = 0$.

Suppose $z = {\beta}_1$. Then, $v(z; z, {\sigma}_1, 0) = 0$, $v(x; z, {\sigma}_1, 0) = +\infty$ for all $x\in[{\alpha}_1, {\beta}_1)$, and $v(x; z, {\sigma}_1, 0) = \ell_1(x)$ for all $x \in [0, {\alpha}_1)$, where $\ell_1$ is the unique positive solution of the singular problem

Observe that, although for this particular problem the uniqueness of $\ell_1$ can be inferred from the phase-portrait 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].

Suppose $z\in ({\beta}_1, {\beta}_2)$. Then, $v(x; z, {\sigma}_1, 0) = \ell_1(x)$ for all $x\in [0, {\alpha}_1)$, $v(x; z, {\sigma}_1, 0) = +\infty$ for every $x \in [{\alpha}_1, {\beta}_1]$, and $v(x; z, {\sigma}_1, 0) = m_1(x)$ for each $x\in ({\beta}_1, z]$, where $m_1(x)$ is the unique positive solution of the singular problem

Suppose $z = {\beta}_2$. Then, $v(x; z, {\sigma}_1, 0) = \ell_1(x)$ if $x\in [0, {\alpha}_1)$, $v(x; z, {\sigma}_1, 0) = +\infty$ if $x \in [{\alpha}_1, {\beta}_1]$, $v(x; z, {\sigma}_1, 0) = \ell_2(x)$ if $x\in ({\beta}_1, {\alpha}_2)$, where $\ell_2(x)$ is the unique positive solution of the singular problem

$v(x; z, {\sigma}_1, 0) = +\infty$ for all $x\in [{\alpha}_2, {\beta}_2)$, and $v(z; z, {\sigma}_1, 0) = 0$.

More generally, suppose $\kappa\geq 3$, $i\in\{2, ..., \kappa-1\}$, and $z\in ({\beta}_{i}, {\beta}_{i+1})$. Then,

where, for every $j\in\{2, ..., \kappa\}$, $\ell_j$ is the unique positive solution of

and, for every $i\in\{2, ..., \kappa\}$, $m_i(x)$ denotes the unique positive solution of

Now, suppose that $z = {\beta}_\kappa$. Then,

Finally, for every $z\in ({\beta}_\kappa, 1]$,

Similarly, for every $z\in ({\alpha}_\kappa, 1)$, $w(x; z, {\sigma}_1, 0)$ is the unique negative solution of (2.3) at ${\lambda} = {\sigma}_1$ and ${\varepsilon} = 0$. Moreover, when $z = {\alpha}_\kappa$, we have that $w(z; z, {\sigma}_1, 0) = 0$, $w(x; z, {\sigma}_1, 0) = -\infty$ for all $x\in ({\alpha}_\kappa, {\beta}_\kappa]$, and $w(x; z, {\sigma}_1, 0) = n_1(x)$ for all $x\in ({\beta}_\kappa, 1]$, where $n_1$ denotes the unique negative solution of the singular problem

And so on... Naturally, the following definition is consistent

and it turns out that ${\varphi}(z, {\sigma}_1, 0) = \lim_{{\lambda}{\uparrow}{\sigma}_1}{\varphi}(z, {\lambda}, 0)$. Moreover, the next result holds; it is [1,Proposition 4.1].

Theorem 2.1. Suppose that $a = a_{\varepsilon}$ is given by $(1.3)$ with ${\varepsilon} = 0$. Then, for every $j\in\{1, ..., \kappa-1\}$ and $z\in [{\beta}_j, {\alpha}_{j+1}]$, we have that

Moreover, for every $j\in\{1, ..., \kappa\}$, there exists $\zeta_j\in({\alpha}_j, {\beta}_j)$ such that

When ${\kappa}\geq 3$, by using the symmetry of $a_{\varepsilon}$ and a phase portrait analysis, it can be shown that, in (2.7), $\zeta_j = \frac{{\alpha}_j+{\beta}_j}{2}$ for every $j\in\{2, ..., \kappa-1\}$. As a consequence of the previous theorem, since ${\varphi}(z, {\lambda}, 0)$ varies continuously with respect to ${\lambda}$ for ${\lambda}\in ({\lambda}_2, {\sigma}_1]$, it is apparent that the next multiplicity result holds; it is Corollary 4.1 of López-Gómez and Rabinowitz ^[1].

Theorem 2.2. Suppose that $a = a_{\varepsilon}$ is given by $(1.3)$ with ${\varepsilon} = 0$. Then, there exists $\eta > 0$ such that, for each ${\lambda}\in [{\sigma}_1-\eta, {\sigma}_1)$, the problem $(1.1)$ possesses, at least, $2\kappa-1$ solutions, $({\lambda}, u_j)$, $1\leq j\leq 2\kappa-1$, having one interior node and such that $u'(0) > 0$.

By the symmetry of the problem, under the assumptions of Theorem 2.2, (1.1) has another $2\kappa-1$ solutions, $({\lambda}, U_j)$, with $U_j'(0) < 0$, $1\leq j\leq 2\kappa-1$. Actually, $U_j = -u_j$, $1\leq j\leq 2\kappa-1$. Therefore, altogether (1.1) posseses, at least, $4\kappa-2$ solutions having one interior node for this range of ${\lambda}$'s.

3.   A multiplicity result for ${\lambda} > \sigma_1$ and sufficiently small ${\varepsilon} > 0$

Naturally, as soon as ${\varepsilon} > 0$, the problem (1.1) with $a = a_{\varepsilon}$ given by (1.3) is non-degenerate, because $a_{\varepsilon}(x) > 0$ for all $x\in [0, 1]$. Thus, the functions $v(x; z, {\lambda}, {\varepsilon})$ and $w(x; z, {\lambda}, {\varepsilon})$ are well defined for every ${\lambda} > {\lambda}_2$ and ${\varepsilon} \in (0, 1]$ if (2.11) holds. Moreover, by uniqueness, they vary continuously with respect to ${\lambda} > {\lambda}_2$, ${\varepsilon} > 0$ and $z\in (0, 1)$. In particular, much like in the degenerate case when ${\varepsilon} = 0$ and ${\lambda}\in ({\lambda}_2, {\sigma}_1]$, we can define ${\varphi}(x; z, {\lambda}, {\varepsilon})$ through (2.1) for these values of ${\lambda}$, ${\varepsilon}$ and $z$. By analyzing the behavior of $v(x; z, {\lambda}, {\varepsilon})$ and $w(x; z, {\lambda}, {\varepsilon})$ as ${\varepsilon}{\downarrow} 0$, we are able to obtain the result holds. It is the first existing multiplicity result for 1-node solutions in the classical non-degenerate setting. Although it resembles Theorem 2.2, in the sense that for some range of ${\lambda}$'s the number of obtained solutions is the same, it is relevant to keep in mind that the nature of this result is completely different, since it deals with the case ${\lambda} > \sigma_1$ and the non-degenerate case ${\varepsilon} > 0$.

Theorem 3.1. Suppose that $a = a_{\varepsilon}$ satisfies $(1.3)$. Then:

(a) for every ${\lambda} > {\sigma}_1$, there exists ${\varepsilon}_0 = {\varepsilon}_0({\lambda}) > 0$ such that, for each ${\varepsilon} \in (0, {\varepsilon}_0]$, the problem $(1.1)$ possesses, at least, $2\kappa-2$ solutions with one interior node;

(b) for every ${\lambda}\in ({\sigma}_1, {\sigma}_2)$, there exists ${\varepsilon}_1 = {\varepsilon}_1({\lambda}) > 0$ such that, for each ${\varepsilon} \in (0, {\varepsilon}_1]$, the problem $(1.1)$ has, at least, $4\kappa-2$ solutions with one interior node.

Remark 3.2. If ${\varepsilon} > 0$ in $\kappa-1$ intervals $({\alpha}_j, {\beta}_j)$, while ${\varepsilon} = 0$ in the remaining one, then, thanks to Theorem 2.3 of López-Gómez and Rabinowitz ^[16], problem (1.1) admits exactly two nodal solutions for each ${\lambda}\in [{\sigma}_1, {\sigma}_2)$. Thus, increasing ${\varepsilon}$ in one single well has dramatic consequences, as, due to Theorem 3.1, (1.1) admits $4\kappa-2$ solutions with one interior node, regardless the value of $\kappa\geq 2$.

In the rest of this section we will provide a proof of Theorem 3.1. We will first show that, for every ${\lambda} > {\sigma}_1$ and sufficiently small ${\varepsilon} > 0$, the function ${\varphi}(z, {\lambda}, {\varepsilon})$ defined by (2.1), with ${\lambda} > {\sigma}_1$, possesses $\kappa-1$ zeroes. Then, we will prove that for each ${\lambda}\in({\sigma}_1, {\sigma}_2)$ ${\varphi}(z, {\lambda}, {\varepsilon})$ possesses $\kappa$ additional zeros. We will derive these properties from a sharp asymptotic analysis of $v(x; z, {\lambda}, {\varepsilon})$ and $w(x; z, {\lambda}, {\varepsilon})$ as ${\varepsilon}{\downarrow} 0$.

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\geq {\beta}_j$, $j\geq 1$, the function $v(\cdot; z, {\lambda}, {\varepsilon})$ blows up in $({\alpha}_i, {\beta}_i)$ for all $i\in\{1, ..., j\}$ as ${\varepsilon}{\downarrow} 0$. It extends some previous findings of ^[24] in a different context.

Note that the map ${\varepsilon}\mapsto v(\cdot; z, {\lambda}, {\varepsilon})$ is decreasing, because, due to [23,Theorem 7.10], the smaller ${\varepsilon}$, the bigger the solution of $(2.2)$. Consequently, its point-wise limit as ${\varepsilon}{\downarrow} 0$ is well defined.

Lemma 3.3. Suppose that ${\lambda} > {\sigma}_1$ and $z\geq {\beta}_j$ for some $j\geq 1$. Then, for every (sufficiently small) $\eta > 0$ and $i \in\{1, ..., j\}$, the solution $v$ of $(2.2)$ satisfies

Proof. Let $i\geq 1$ be an integer such that $i\leq j$ and set $V: = \sqrt{{\varepsilon}}v$. Then, in the interval $({\alpha}_i, {\beta}_i)$, $V$ solves the differential equation $-V'' = {\lambda} V-V^3$. Thus, since $V({\alpha}_i) > 0$ and $V({\beta}_i)\geq 0$, $V$ provides us with a positive supersolution of the homogeneous boundary value problem

Since ${\lambda} > {\sigma}_1$, this problem admits a unique positive solution, denoted by $V_0$. By [23,Th. 7.10], it follows that $V_0\leq V = \sqrt{{\varepsilon}}v$. Finally, (3.1) is a direct consequence of the estimate $v\geq \frac{V_0}{\sqrt{{\varepsilon}}}$ by letting ${\varepsilon}{\downarrow} 0$.

With some additional effort, and using the specific structure of the weight $a_{\varepsilon}$ given by (1.3), we are able to obtain the following result.

Proposition 3.4. Suppose that ${\lambda} > {\sigma}_1$ and $z > {\beta}_j$ for some $j\geq 1$. Then, for every $i \in\{1, ..., j\}$,

Proof. Thanks to Lemma 3.3, it suffices to prove

For convenience, we set $v_{\varepsilon}: = v(\cdot; z, {\lambda}, {\varepsilon})$, and divide the proof into several steps.

Step 1: proof of (3.3) for ${\alpha_1}$. Since $a\equiv 1$ in $[0, {\alpha}_1]$, it is easily seen that, for every $x\in (0, {\alpha}_1)$ and ${\varepsilon}\in (0, 1]$,

Based on [5,Th. 3.4], the unique positive solution, $\ell_1$, of the singular problem

satisfies $\ell_1 = \lim_{M{\uparrow}\infty}L_M$, where $L_M$ is the unique positive solution of

The uniqueness of $\ell_1$, as well as of the remaining large positive solutions of this paper, is a direct consequence of the theory of López-Gómez and Maire ^[21] and López-Gómez, Maire and Véron ^[22]. Then, thanks to [23,Th. 7.10], it is easily seen that, for sufficiently large $M > 0$, $v_{\varepsilon}\leq L_M \leq \ell_1$ in $[0, {\alpha}_1)$. Hence, $v_0: = \lim_{{\varepsilon}{\downarrow} 0}v_{\varepsilon} \leq \ell_1$ (observe that this point-wise limit is well defined by the monotonicity of $v_{\varepsilon}$ with respect to ${\varepsilon}$). Moreover, with slightly more effort, a standard phase portrait analysis also shows that $v_{\varepsilon}'(0) < \ell_1'(0)$ if ${\varepsilon}\in (0, 1]$, as otherwise $v_{\varepsilon}$ would not be defined in the whole interval $[0, {\alpha}_1]$. Thus, by (3.5),

for all ${\varepsilon}\in (0, 1]$. In order to prove that $\lim_{{\varepsilon}{\downarrow} 0} v_{\varepsilon}({\alpha}_1) = +\infty$, we can argue by contradiction assuming that there exists a constant $C > 0$ such that $v_{\varepsilon}({\alpha}_1)\leq C$ for all ${\varepsilon}\in (0, 1]$. Then, by (3.5), there also exists a positive constant $\tilde C$ such that $|v_{\varepsilon}'({\alpha}_1)|\leq \tilde C$ for all ${\varepsilon}\in (0, 1]$, and, setting $V_{\varepsilon}: = \sqrt{{\varepsilon}}v_{\varepsilon}$, it is apparent that

Since $-V_{\varepsilon}'' = {\lambda} V_{\varepsilon}-V_{\varepsilon}^3$ in $[{\alpha}_1, {\beta}_1]$, it follows from (3.6) that $V_{\varepsilon}\to 0$ in $[{\alpha}_1, {\beta}_1]$ as ${\varepsilon}{\downarrow} 0$, which contradicts $V_{\varepsilon}\geq V_0 > 0$, where $V_0$ is the unique positive solution of (3.2). Therefore, $v_{\varepsilon}({\alpha}_1){\uparrow} +\infty$ as ${\varepsilon} {\downarrow} 0$. Note that, since $L_{v_{\varepsilon}({\alpha}_1)} = v_{\varepsilon}$ in $[0, {\alpha}_1]$, we can infer from [5,Th. 3.4] that

Step 2: proof of $(3.3)$ for ${\beta_1}$. Since $a_{\varepsilon}\equiv {\varepsilon}$ in $[{\alpha}_1, {\beta}_1]$, owing to (3.5), it becomes apparent that, for every ${\varepsilon}\in (0, 1]$,

Thus, by Step 1, we find that

To show that $\lim_{{\varepsilon}{\downarrow} 0}v_{\varepsilon}({\beta}_1) = +\infty$, we will argue by contradiction assuming that $\{v_{\varepsilon}({\beta}_1)\}_{{\varepsilon}\in (0, 1]}$, is bounded from above. We first show that, for ${\varepsilon}\sim 0$, $v_{\varepsilon}$ is decreasing in a left neighborhood of ${\beta}_1$. This will imply, together with (3.7), that

Indeed, in the contrary case, by Lemma 3.3 and our contradiction hypothesis, for every ${\varepsilon}\sim 0$ there would exist $y_{\varepsilon} < {\beta}_1$, $y_{\varepsilon}\sim{\beta}_1$, such that where $v_{\varepsilon}(y_{\varepsilon}) < v_{\varepsilon}({\beta}_1)$ and $y_{\varepsilon}$ is a local minimum point. Thus,

which implies $0\geq \lambda-{\varepsilon} v_{\varepsilon}^2(y_{\varepsilon})$. Thus, letting ${\varepsilon}{\downarrow} 0$, we find that ${\lambda}\leq 0$, which contradicts ${\lambda} > {\sigma}_1$ and concludes the proof of (3.8).

Let $\eta > 0$ be such that ${\beta}_1+\eta < \min\{z, {\alpha}_2\}$. By the monotonicity of $v_{\varepsilon}$ with respect to ${\varepsilon}$, either $\{v_{\varepsilon}\}_{{\varepsilon}\in (0, 1]}$ is bounded in $[{\beta}_1, {\beta}_1+\eta]$, or

Suppose the first option occurs. Then, it follows from the differential equation that $v_{\varepsilon}''$ also is bounded in $[{\beta}_1, {\beta}_1+\eta]$ uniformly in ${\varepsilon}\in (0, 1]$. Thus, since

and $\int_{{\beta}_1}^x v_{\varepsilon}''(s)\, ds$ is bounded, we can infer from (3.8) that $\lim_{{\varepsilon}{\downarrow} 0}v_{\varepsilon}'(x) = -\infty$ uniformly in $x\in [{\beta}_1, {\beta}_1+\eta]$. Consequently, for every integer $n\geq 1$, there exists ${\varepsilon}_0 = {\varepsilon}_0(n) > 0$ such that $v_{\varepsilon}'\leq -n$ in $[{\beta}_1, {\beta}_1+\eta]$ if ${\varepsilon} < {\varepsilon}_0$. Hence,

for all ${\varepsilon} < {\varepsilon}_0$. Therefore, since $v_{\varepsilon}({\beta}_1)$ is bounded and $n$ can be is arbitrarily large, we find that $v_{\varepsilon}({\beta}_1+\eta) < 0$ for sufficiently small ${\varepsilon} > 0$, which is impossible, because $v_{\varepsilon}$ is positive everywhere. This contradiction shows that (3.9) holds.

For every ${\varepsilon}\in (0, 1]$, let $x_{\varepsilon}\in [{\beta}_1, {\beta}_1+\eta]$ be such that $v_{\varepsilon}(x_{\varepsilon}) = \|v_{\varepsilon}\|_{ \mathcal {C}[\beta_1, \beta_1+\eta]}$. Since $v_{\varepsilon}({\beta}_1)$ is bounded, (3.12) implies that $x_{\varepsilon}\in ({\beta}_1, {\beta}_1+\eta]$ for sufficiently small ${\varepsilon} > 0$. Suppose that $x_{\varepsilon}\in ({\beta}_1, {\beta}_1+\eta)$. Then, $v_{\varepsilon}'(x_{\varepsilon}) = 0$ and $v_{\varepsilon}''(x_{\varepsilon})\leq 0$. Thus,

which implies $v_{\varepsilon}(x_{\varepsilon})\leq \sqrt{{\lambda}}$ and contradicts (3.9). Therefore, $x_{\varepsilon} = {\beta}_1+\eta$ for sufficiently small ${\varepsilon} > 0$. So, according to (3.9), we find that $\lim_{{\varepsilon}{\downarrow} 0}v_{\varepsilon}({\beta}_1+\eta) = +\infty$, which contradicts the estimate $v_{\varepsilon}\leq \ell_2$, where $\ell_2$ denotes the unique positive solution of the singular problem

This contradiction shows that $\lim_{{\varepsilon}{\downarrow} 0}v_{\varepsilon}({\beta}_1) = +\infty$.

Step 3: proof of $(3.3)$ in the remaining cases. Now, we suppose that $z > {\beta}_2$ and set $v_0({\alpha}_2): = \lim_{{\varepsilon}{\downarrow} 0}v_{\varepsilon}({\alpha}_2)$. This limit (which is either a positive real number, or $+\infty$) exists by the monotonicity of $v_{\varepsilon}$ with respect to ${\varepsilon}\in (0, 1]$. Moreover, by [23,Th. 7.10], $v_{\varepsilon}\leq \ell_2$ in $[{\beta}_1, {\alpha}_2]$. Thus, by a standard compactness argument, the point-wise limit $v_0: = \lim_{{\varepsilon}{\downarrow} 0} v_{\varepsilon}$ in $[{\beta}_1, {\alpha}_2]$ provides us with the unique positive solution of the singular problem

Actually, $v_{\varepsilon}\to v_0$ in $\mathcal {C}^1[{\beta}_1+{\delta}, {\alpha}_2-{\delta}]$ as ${\varepsilon}{\downarrow} 0$ for all ${\delta} > 0$, ${\delta}\sim 0$. Thus, since for every $x\in ({\beta}_1, {\alpha}_2)$ and ${\varepsilon} > 0$, we have that

it becomes apparent that the right hand side of (3.10) remains bounded as ${\varepsilon}{\downarrow} 0$. Therefore, by adapting the argument of proof of Step 1, it readily follows that $v_0({\alpha}_2) = +\infty$. As a byproduct, $v_0 = \ell_2$ in $({\beta}_1, {\alpha}_2)$. Furthermore, since

the proof of Step 2 can be as well adapted to conclude that

This procedure can be repeated for the other ${\alpha}_i$'s and ${\beta}_i$'s, and ends the proof.

Remark 3.5. From the proof of Proposition 3.4, it is apparent that, as soon as $z > {\beta}_{j}$ for some $j\geq 2$, $\lim_{{\varepsilon}{\downarrow} 0}v_{\varepsilon} = \ell_j$ in $[{\beta}_{j-1}, {\alpha}_j]$, where $\ell_j$ is the unique positive solution of

Naturally, the negative solution $w_{\varepsilon}(x): = w(x; z, {\lambda}, {\varepsilon})$ satisfies similar properties as $v_{\varepsilon}$ when ${\varepsilon} {\downarrow} 0$. Thus, to obtain the proof of Theorem 3.1 from Proposition 3.4 and Remark 3.5, one can argue as follows.

Proof of Theorem 3.1. (a) Suppose that $z \in ({\beta}_{j-1}, {\alpha}_j)$ for some $j\in\{2, \dots, {\kappa}\}$. Then, based on the generalized Keller–Osserman condition of ^[27], the a priori bounds of Keller ^[25] and Osserman ^[26], and the point-wise monotonicity of $v_{\varepsilon}(x): = v(x; z, {\lambda}, {\varepsilon})$ as ${\varepsilon}{\downarrow} 0$, it becomes apparent that, in the interval $[{\beta}_{j-1}, z]$, $v_{\varepsilon}$ approximates, as ${\varepsilon}{\downarrow} 0$, the unique positive solution of

denoted by $p_z$. Similarly, in $[z, {\alpha}_j]$, $w_{\varepsilon}$ approaches, as ${\varepsilon}{\downarrow} 0$, the unique negative solution of

denoted by $n_z$. Since a careful phase portrait analysis reveals that

by continuous dependence, it becomes apparent that, for sufficiently small ${\varepsilon} > 0$,

Therefore, for every $j\in\{2, \dots, {\kappa}\}$, ${\varphi}(z({\varepsilon}), {\lambda}, {\varepsilon}) = 0$ for some $z({\varepsilon})\in ({\beta}_{j-1}, {\alpha}_j)$, which provides us with at least $2(\kappa-1)$ different solutions of (1.1) with one interior node.

(b) When ${\lambda}\in({\sigma}_1, {\sigma}_2)$, to construct the remaining $2\kappa$ solutions with one node, we work in $({\alpha}_i, {\beta}_i)$ for any given $i\in\{1, ..., \kappa\}$. Once $\lambda$ has been fixed, we define

Since ${\sigma}_1 < {\lambda} < {\sigma}_2$, the following relations hold

Thus, for every $z\in (z_{-, i}({\lambda}), z_{+, i}({\lambda}))$, we have that

Consequently, by [5,Th. 4.1], the functions $v_{\varepsilon}: = v(\cdot; z, {\lambda}, {\varepsilon})$ and $w_{\varepsilon}: = w(\cdot; z, {\lambda}, {\varepsilon})$ are well defined in $({\alpha}_i, z]$ and $[z, {\beta}_i)$, respectively, for all ${\varepsilon} \in [0, 1]$, as well as

We start by considering ${\varepsilon} = 0$. As $z < \tfrac{{\alpha}_i+{\beta}_i}{2}$ for $z\sim z_{i}^{-}({\lambda})$,

Thus, [28,Th. 2.4] gives

With some additional effort, one can prove that $\lim_{z{\downarrow} z_{i}^{-}({\lambda})}w_{0}'(z) = -\infty$, because $w_{0}(z) = 0$. On the other hand, the inequality in (3.24) guarantees that $v_0$ converges to some finite profile in $({\alpha}_i, z]$ as $z{\downarrow} z_{i}^{-}({\lambda})$; thus $\lim_{z{\downarrow} z_{i}^{-}({\lambda})}v_{0}'(z)$ is finite. Hence, we can conclude that ${\varphi}(z, {\lambda}, 0) > 0$ provided $z > z_{i}^{-}({\lambda})$ is sufficiently close to $z_{i}^{-}({\lambda})$. By continuity, ${\varphi}(z, {\lambda}, {\varepsilon}) > 0$ for sufficiently small ${\varepsilon} > 0$.

Similarly, one can prove that ${\varphi}(z, {\lambda}, {\varepsilon}) < 0$ for sufficiently small ${\varepsilon} > 0$ provided $z < z_{i}^{+}({\lambda})$ is sufficiently close to $z_{i}^{+}({\lambda})$. Therefore, ${\varphi}$ admits $\kappa$ zeroes which are different from the ones obtained in part (a), since they lie in $({\alpha}_i, {\beta}_i)$, $i\in\{1, \dots, {\kappa}\}$. Equivalently, (1.1) admits, at least, $2\kappa$ solutions with one node for sufficiently small ${\varepsilon} > 0$, besides the $2\kappa-2$ solutions constructed in Part (a).

4.   An example with $\kappa=2$

In this section, we provide some numerical results related to the case (1.3) with $\kappa = 2$,

and ${\varepsilon}\in [0, 1)$. Then,

Figure 5 shows the global bifurcation diagram of 1-node solutions, $({\lambda}, u)$, with $u'(0) > 0$, for ${\varepsilon} = 0$. As remarked above, once the sign of the derivative at $x = 0$ has been chosen, every solution $({\lambda}, u)$ is uniquely identified, apart from the value of ${\lambda}$ of course, by the position of its interior zero $z$. Thus, in all the bifurcation diagrams of this paper, we represent the values of $z$ in abscissas versus the values of ${\lambda}$ in ordinates.

In agreement with Theorem 2.2, (1.1) has three 1-node solutions with $u'(0) > 0$ for each ${\lambda} < {\sigma}_1$ sufficiently close to ${\sigma}_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\sim 0.35$ and $z\sim 0.65$, respectively, if ${\lambda}\sim{\sigma}_1$. Observe that those values are the midpoints of the vanishing intervals of $a_{0}(x)$, $[0.3, 0.4]$ and $[0.6, 0.7]$, respectively. These branches bifurcate from the branch of 1-node solutions with $z = 0.5$ at ${\lambda}_b\sim 104.8215$.

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 $({\lambda}, u) = ({\lambda}, 0)$ at ${\lambda} = {\lambda}_2$. Thanks to the exchange stability principle ^[29], it consists of unstable solutions until ${\lambda}$ reaches ${\lambda}_b$. At such a point, the solutions with $z = 0.5$ become stable for any further value of ${\lambda}$. Once again, by the exchange stability principle, the curves that bifurcate from the branch of solutions with $z = 0.5$ consist of unstable solutions until ${\lambda}$ reaches ${\sigma}_1$ where they become metasolutions, much like $v$ and $w$ introduced in (2.2) and (2.3), respectively, as ${\lambda}{\uparrow} {\sigma}_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 ${\lambda}$. In the first picture we have plotted the nodal solutions with $z = 0.5$. It is apparent that, as ${\lambda}$ increases, they become much larger in

than in the complement of $a^{-1}_{0}(0)$. Actually, as already discussed in Section 2, these solutions blow up, as ${\lambda}{\uparrow} {\sigma}_1$, to $+\infty$ in $[0.3, 0.4]$, and to $-\infty$ in $[0.6, 0.7]$, while they stabilize to the appropriate large solutions in the complement of $a^{-1}_{0}(0)$, in agreement with [1,Theorem 2.2]. In particular, in the interval $[0.4, 0.6]$ they converge, as ${\lambda}{\uparrow} {\sigma}_1$, to the unique 1-node solution of the singular problem

The central plot of Figure 5 shows a series of solutions on the left bifurcated branch of Figure 4. The node of these solutions approaches $0.35$, the middle point of the interval $[0.3, 0.4]$, as ${\lambda}{\uparrow} {\sigma}_1$. These solutions blow-up to $-\infty$ in $[0.6, 0.7]$ as ${\lambda}{\uparrow} {\sigma}_1$, while in the complement they stabilize to certain large solutions, confirming what is established by [1,Theorem 2.2]. In particular, in the interval $[0, 0.6]$ they converge, as ${\lambda}{\uparrow} {\sigma}_1$, to the unique large solution with one interior node at $0.35$ of the singular problem

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 $+\infty$ in $[0.3, 0.4]$ as ${\lambda}{\uparrow}{\sigma}_1$.

As expected from Theorem 3.1(a), for sufficiently small ${\varepsilon} > 0$, some of the previous nodal solutions are defined for all values of ${\lambda} > {\sigma}_1$. Surprisingly, our numerical computations show that this is true also for ${\varepsilon}$ close to 1. For example, Figure 7 shows the corresponding global bifurcation diagram for ${\varepsilon} = 0.9999$.

Essentially, the solutions behave in the same way as in the special case when ${\varepsilon} = 0$, except for the relevant fact that they are defined for all values of ${\lambda} > {\sigma}_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 ${\lambda}$ increases, and the third one superimposes the plots of a series of solutions on the right secondary branch, whose nodes approach $0.65$ as ${\lambda}$ 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 ${\varepsilon} > 0$, the existence of one solution for every ${\lambda} > {\sigma}_1$ and of three if ${\lambda}\in({\sigma}_1, {\sigma}_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 ${\lambda}$ and ${\varepsilon}$ 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 ${\varepsilon} = 0.99999999$, though in such a case the bifurcation point from the central branch grew up to ${\lambda}_b = 435.0634$. This causes a great perplexity to us, as problem (1.1) with $a\equiv 1$ has a unique solution with one interior node at $z = 0.5$, whereas the nodes of the solutions along the secondary branches for ${\varepsilon} = 0.99999999$ still approach $0.35$ and $0.65$ as ${\lambda}$ increases for really huge intervals of ${\lambda}$. The extremely challenging problem of ascertaining how these nodal solutions disappear as ${\varepsilon}{\uparrow} 1$, and only one solution remains for ${\varepsilon} = 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 ${\lambda}$, and that such a loop shrinks to a single point up to disappear as ${\varepsilon} {\uparrow} 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 ${\lambda}_b({\varepsilon})$ converges to $\infty$ as ${\varepsilon}{\uparrow} 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 ${\lambda}$, we are working close to the computational limits of the computer.

5.   An example with $\kappa=3$.

In this section we present the results of our numerical computation related to a weight $a_{{\varepsilon}}$ as in (1.3) with $\kappa = 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 ${\varepsilon} = 0$. In agreement with Theorem 2.2, problem (1.1) has 5 nodal solutions $({\lambda}, u)$ with $u'(0) > 0$ for ${\lambda} < {\sigma}_1$ sufficiently close to ${\sigma}_1$. Nevertheless, our computations show that there is some intermediate range of values of ${\lambda}$ 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 $\mathcal{G}_2^+$, plus two additional components: one on the left of $\mathcal{G}_2^+$, say $\mathcal{H}_{2, \ell}^+$, and another one on its right, $\mathcal{H}_{2, r}^+$, which are isolated with respect to the trivial solution $({\lambda}, u) = ({\lambda}, 0)$.

Moreover, observe that $\mathcal {G}_2^+$ contains a closed loop of unstable solutions. It surrounds the solutions with $z = 0.5$ and bifurcates from such a curve at ${\lambda}_{b, 1} = 66.4542$ and ${\lambda}_{b, 2} = 207.4892$. The two lateral folding type components have their turning points at

Thus, for every ${\lambda}\in ({\lambda}_t, {\lambda}_{b, 2})$ problem (1.1) has, at least, $14$ solutions, $({\lambda}, u)$: $7$ among them satisfy $u'(0) > 0$, and their opposite, with $u'(0) < 0$. The solutions with $z = 0.5$ are unstable for ${\lambda}\in ({\lambda}_2, {\lambda}_{b, 1})$, stable for ${\lambda}\in ({\lambda}_{b, 1}, {\lambda}_{b, 2})$ and unstable for ${\lambda}\in ({\lambda}_{2, b}, {\sigma}_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 $\mathcal {G}_2^+$. The solutions on the first picture were computed for ${\lambda}$ close to ${\lambda}_2 = (2\pi)^2$, while the solutions on the second picture where computed for a range of ${\lambda}$'s closer to ${\sigma}_1$. These solutions blow up to $+\infty$ in $[0.2, 0.3]$, and to $-\infty$ in $[0.7, 0.8]$, as ${\lambda}{\uparrow} {\sigma}_1$, whereas they stabilize on the complement of these two intervals to the appropriate large solutions, either positive, negative, or nodal.

Figure 11 shows the plots of a series of solutions with one interior node along the closed loop of $\mathcal {G}_2^+$. Although the nodes of these solutions perturb from $z = 0.5$ at the bifurcation values ${\lambda}_{b, 1}$, ${\lambda}_{b, 2}$, and separate away from $0.5$ a substantial amount, reaching approximatively $0.4$ in the first picture and $0.6$ in the second one, observe that they do not reach the remaining two intervals of $a^{-1}_{0}(0)$.

Figure 12 shows the plots of a series of solutions on each of the two half-branches of the component $\mathcal {H}_{2, \ell}^+$ for a series of values of ${\lambda}$ up to 400. The solutions on the first picture, where $z\sim 0.25$, are unstable, while those on the second picture are stable and satisfy $z\sim 0.375$ (see Figure 9). Note that $0.375$ is the midpoint of $[0.3, 0.45]$, in agreement with Theorems 2.1 and 2.2. Moreover, our numerical simulations show that also when the node $z$ lies in $[0.2, 0.3]$, it converges to the midpoint of such interval. This actually goes beyond the scope of Theorems 2.1 and 2.2. The solutions of the first picture must blow-up to $-\infty$ in $[0.45, 0.55]\cup [0.7, 0.8]$ as ${\lambda}{\uparrow} {\sigma}_1$, while in the complement of these two intervals they stabilize to the appropriate large solutions. The solutions on the second picture, instead, must blow-up to $+\infty$ in $[0.2, 0.3]$ and to $-\infty$ in $[0.45, 0.55]\cup[0.7, 0.8]$, while they stabilize on $[0, 1]\setminus a^{-1}_{0}(0)$ to the appropriate large solutions, either positive, negative, or nodal, as ${\lambda}{\uparrow} {\sigma}_1$.

Finally, Figure 13 shows the plots of a series of nodal solutions on the component $\mathcal {H}_{2, r}^+$ for the same range of values of ${\lambda}$ used in Figure 12. The profiles are similar to those computed along $\mathcal {H}_{2, \ell}^+$, except for the fact that now the nodes approach the values $0.625$ and $0.75$, respectively, as ${\lambda}{\uparrow} {\sigma}_1$. In particular, these solutions blow up to $+\infty$ in $[0.2, 0.3]\cup [0.45, 0.55]$ as ${\lambda}{\uparrow} {\sigma}_1$. Moreover, the stable solutions, those with $z\sim 0.625$, also blow up to $-\infty$ in $[0, 7, 0.8]$ as ${\lambda}{\uparrow} {\sigma}_1$.

This precise point-wise behavior of the nodal solutions as ${\lambda}{\uparrow} {\sigma}_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 $\mathcal {G}_2^+$ (or, equivalently, $\mathcal {C}_2^+$), we numerically solved the bifurcation problem from $u = 0$ at ${\lambda} = {\lambda}_2$ and then continued the bifurcated curve, detecting and marking the two bifurcation points along it. We stopped when we reached a value of ${\lambda}$ sufficiently close to ${\sigma}_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.

As far as the global folding type components are concerned, as they do not bifurcate from $u = 0$ but from certain metasolutions at ${\lambda} = {\sigma}_1$, their existence must be inferred from Theorem 2.2. Indeed, thanks to Theorem 2.2, (1.1) has five solutions for ${\lambda} < {\sigma}_1$, ${\lambda}\sim {\sigma}_1$, and, thanks to ^[1,7,16], their sharp point-wise behavior as ${\lambda}{\uparrow} {\sigma}_1$ is well known (see Figures 10, 12 and 13). The knowledge of such limiting behavior has been crucial to compute the components $\mathcal {H}_{2, \ell}^+$ and $\mathcal {H}_{2, r}^+$, as described in Section 7.

Once we completed this analysis, we considered the non-degenerate problem (1.1) with the choice (5.1) and ${\varepsilon}\in (0, 1)$. Figure 14 shows the corresponding global bifurcation diagrams for ${\varepsilon} = 0.1$, ${\varepsilon} = 0.7$, ${\varepsilon} = 0.8$ and ${\varepsilon} = 0.9999$, respectively.

In the global bifurcation diagram for ${\varepsilon} = 0.1$, the loop of solutions on $\mathcal {G}_2^+$ bifurcates from the solution with $z = 0.5$ at ${\lambda}_{b, 1} = 69.5268$ and ${\lambda}_{b, 2} = 193.3319$, while the turning points of the components $\mathcal {H}_{2, \ell}^+$ and $\mathcal {H}_{2, r}^+$ occur at ${\lambda}_t = 173.5823$. In the global bifurcation diagram for ${\varepsilon} = 0.7$, these values changed to ${\lambda}_{b, 1} = 110.7253$, ${\lambda}_{b, 2} = 160.4316$, and ${\lambda}_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 ${\varepsilon}$ until some critical value ${\varepsilon}_c\in (0.77, 0.78)$ where it shrinks to a single point before disappearing for all further values of ${\varepsilon}$ 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 ${\varepsilon} = 0.8$, ${\lambda}_t = 230.2345$. Finally, ${\lambda}_t = 1.1068\times 10^3$ when ${\varepsilon} = 0.9999$. Thus, it is apparent that, although ${\lambda}_t$ increases as the secondary bifurcation parameter ${\varepsilon}$ 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 ${\varepsilon}$ for which we have computed them, going much beyond the local perturbation result given by Theorem 3.1.

6.   An example with $\kappa=4$

In this section we will discuss the model (1.1) for the choice (1.3) with $\kappa = 4$, ${\varepsilon}\in[0, 1)$ and (1.5). According to Theorems 2.2 and 3.1 the associated problems admit, at least, $14$ solutions with one interior node, for ${\lambda}$ in a left neighborhood of ${\sigma}_1$ when ${\varepsilon} = 0$ and for ${\lambda}\in({\sigma}_1, {\sigma}_2)$ and ${\varepsilon} > 0$ sufficiently small, respectively. Moreover, for ${\lambda}\geq{\sigma}_2$ and ${\varepsilon} > 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 ${\varepsilon} = 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 $({\lambda}, u)$ with one interior node and $u'(0) > 0$. By the symmetry of the problem, the pairs $({\lambda}, -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, $\mathcal {G}_2^+$ (or, equivalently, $\mathcal {C}_2^+$), bifurcating from $u = 0$ at ${\lambda} = {\lambda}_2$, plus two additional supercritical folds, $\mathcal {H}_{2, \ell}^+$ and $\mathcal {H}_{2, r}^+$. The component $\mathcal {G}_2^+$ consists of the \lq"symmetric" solutions with $z = 0.5$ plus a secondary curve of solutions bifurcating supercritically from the primary curve at ${\lambda}_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 $\mathcal {H}_{2, \ell}^+$ and $\mathcal {H}_{2, r}^+$ occur at ${\lambda}_t = 433.9062 > {\lambda}_b$, problem (1.1) has --- limited to solutions with one interior node and $u'(0) > 0$ ---, at least, one solution $({\lambda}, u)$ for each ${\lambda}\in({\lambda}_2, {\lambda}_b]$, three solutions for each ${\lambda}\in ({\lambda}_b, {\lambda}_t)$, five at ${\lambda} = {\lambda}_t$ and seven for each ${\lambda}\in({\lambda}_t, {\sigma}_1)$. The nodes of these solutions, as well as their point-wise behaviors as ${\lambda}{\uparrow} {\sigma}_1$, obey the general patterns described in Section 2. Therefore, the solution blows up to $\pm\infty$ as ${\lambda}{\uparrow}{\sigma}_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 $+\infty$, while it approximates $-\infty$ 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 ${\varepsilon} > 0$, ${\varepsilon}\sim 0$, at least for ${\lambda}\in({\sigma}_1, {\sigma}_2)$. Nevertheless, our numerical experiments show that the same global structure persists for all further values of ${\varepsilon}\in (0, 1)$ for which we have computed it, and the full multiplicity result given by Theorem 3.1(b) also holds for ${\lambda}\geq{\sigma}_2$. Figure 16 shows a series of global bifurcation diagrams computed for ${\varepsilon} = 0.1$, ${\varepsilon} = 0.5$, ${\varepsilon} = 0.9$ and ${\varepsilon} = 0.9999$, confirming that our multiplicity results, of a local nature, are actually valid, at least, for all ${\varepsilon}\in (0, 0.9999]$.

Table 1 provides us with the values of ${\varepsilon}$, ${\lambda}_b$ and ${\lambda}_t$ in each of the global bifurcation diagrams plotted in Figure 15.

As in all cases analyzed in the previous sections, the analysis of the fine behavior of the global bifurcation diagram as ${\varepsilon}{\uparrow} 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 vary with the distribution of the intervals where $a_{{\varepsilon}}(x) = {\varepsilon}$, as well as with the value of the parameter ${\varepsilon}$. In particular, some of our numerical simulations (not presented here) suggest that there are examples where, besides the components $\mathcal {G}_2^+$, $\mathcal {H}_{2, \ell}^+$ and $\mathcal {H}_{2, r}^+$, the model exhibits some additional \lq"isola"\, also containing nodal solutions. A complete description of these more sophisticated aspects will appear elsewhere.

7.   On the numerical experiments

In most of previous works devoted to numerical continuation of reaction-diffusion problems, a pseudo-spectral method combining a trigonometric spectral method with collocation at equidistant points has been usually used to discretize (1.1), especially to compute the small solutions bifurcating from $u = 0$ (see, e.g., Gómez-Reñasco and López-Gómez ^[31], López-Gómez, Eilbeck, Duncan and Molina-Meyer ^[32], López-Gómez and Molina-Meyer ^[33], López-Gómez, Molina-Meyer and Tellini ^[34], López-Gómez, Molina-Meyer and Rabinowitz ^[7], and Fencl and López-Gómez ^[35]).

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 9). 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\leq N$) of the tridiagonal Toeplitz matrix corresponding to the discretization of the operator $-\frac{\mathrm{d}^2}{\mathrm{d}x^2}$ is given by

(see, e.g., [38,Section 2]), and it converges to the exact value ${\lambda}_k = (k\pi)^2$ as $N\to+\infty$. For our computations, we chose $N = 501$, so that

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 $\mathcal{C}_2^+$ (or, equivalently, $\mathcal{G}_2^+$), we started from a value of ${\lambda}$ bigger than, but close to the bifurcation point ${\lambda}_2$ and used the initial iterate $u_0 = \sin(2\pi x)$, since, there, it provides us with a first order approximation of the shape of the solutions bifurcating from $({\lambda}_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 $\mathcal{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 ${\lambda} < \sigma_1$, ${\lambda}\sim \sigma_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 $({\lambda}, 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 ${\lambda}$ close to ${\sigma}_1$ when ${\varepsilon}\sim 0$, as well as for all ${\lambda} > {\lambda}_2 = (2\pi)^2$ when ${\varepsilon}\sim 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 the special choice (1.3) with ${\varepsilon} = 0$, the nodal solutions grew to $+\infty$, or $-\infty$, in some of the intervals $[{\alpha}_j, {\beta}_j]$ as ${\lambda}{\uparrow} {\sigma}_1$, and, actually, the fact that (1.1) has a unique nodal solution $({\lambda}, u)$ with $u'(0) > 0$ for ${\varepsilon} = 1$, while it can admit an arbitrarily large number of nodal solutions with $u'(0) > 0$ for ${\varepsilon}$ really close to $1$, is an amazing (new) phenomenology that could not have been established numerically without the analytical results established in ^[1].

We now detail the computational parameters used in this paper to obtain the bifurcation diagrams. As a stopping criterion for the Newton method, we considered an iteration to be satisfactory when the norm of our augmented system, evaluated at such an iteration, was smaller than $0.0001$. In order to determine bifurcation points, a bisection method with the same tolerance has been implemented, in order to find points where the determinant of the Jacobian of the augmented system vanishes.

Since at some point, solutions became very large in some intervals and, hence, computational times started to increase considerably, an adaptive step in the Keller–Yang algorithm has been implemented. An initial step size equal to $0.3$ has been chosen, allowing then the algorithm to vary the step size between $0.1$ and $30$, depending on the precision required. Close to turning points or when a new branch was initiated, high precision or, equivalently, a very small step size was required, while far from these points a lower precision was acceptable.

An empirical method consisting in detecting the rate of change of the slope of the bifurcation curves with respect to $\lambda$ has provided acceptable results at a very low computational cost. In particular, we used the quantity $u'(0)$ to measure such changes and a threshold value of $0.5$ to decide whether to shorten or to increase the step length. We remark that, in order to avoid computational problems, determining when the bifurcation diagrams become very narrow, for instance at the turning points, is extremely important, as pointed out in ^[34,48].

Acknowledgments

This work has been supported by the Ministry of Science, Technology and Universities of Spain, under Research Grant PGC2018-097104-B-I00, and by the IMI of Complutense University of Madrid.

Conflict of interest

The authors declare there is no conflicts of interest.