GLOBAL BIFURCATION DIAGRAMS OF ONE NODE SOLUTIONS IN A CLASS OF DEGENERATE BOUNDARY VALUE PROBLEMS

. In [12], the structure of the set of possible solutions of a degenerate boundary value problem was studied. For solutions with one interior zero, there were two possibilities for the solution set. In this paper, numerical examples are given showing each of these possibilities can occur.

Then the possible behavior of the set of one node solutions becomes more complicated and is described by Theorem 4.2 of [12], which is stated next. Let P denote the projector of R × C 2 [−L, L] to R, i.e. P(λ, u) = λ.
is a continuous curve parameterized by λ. Moreover, either and where T ± 2 denotes the set of nontrivial solutions, (λ, u), such that ±u (−L) > 0 and u has exactly one interior zero in (−L, L).
By Theorem 1.1, the λ-projection of C 2 might just consist of the interval If this occurs, Theorem 1.1 shows C 2 is unbounded as λ → ( π β−α ) 2 and (1) must possess an additional pair of unbounded solution components each member of which lies in T ± 2 . One of them, D + 2 , consists of solutions (λ, u) with u (0) > 0, and its twin, D − 2 = {(λ, −u) : (λ, u) ∈ D + 2 }. Whether or not each of (i) and (ii) of the theorem can actually occur was left as an open problem in [12]. The main goal of this paper is to provide some analytical and numerical evidence supporting the conjecture that either of the alternatives of Theorem 1.1 and, in particular, C 2 = D 2 are possible.
Towards this end the condition (2) will be weakened to allow a −1 (0) to consist of multiple components. Applying some of the results of [12] that are needed here to the case when a −1 (0) has multiple components requires minor adjustments of some of the proofs of [12], which be carried out in Section 2. Then, in Section 3, an example will be given which provided the motivation for the numerical work of Section 4. Moreover this example together with an additional conjecture allows us to give a heuristic argument showing the set of solutions with one node for the example possesses a component that is unbounded at λ = ( π β−α ) 2 . In Section 4, we will perform a series of numerical experiments on this example to show (A) that in some parameter range, the unbounded component is an isolated curve or sort of isola that is separated from C 2 = D 2 , while (B) in another parameter range, C 2 = D 2 .
2. Some extensions of [12]. The paper [12] mainly treats (1) assuming a −1 (0) consists of a single component, although a few results allow a −1 (0) to be a union of several intervals of equal length. The main goal of this section is to provide extensions of some of the main results of [12] for that new condition on a −1 (0). They will be required in Section 3 to perform the construction of the main example of this paper. Henceforth instead of (2), we will assume that there is a constant, h > 0, an integer n ≥ 2 and n disjoint subintervals Without lost of generality, we can assume that β j < α j+1 for all 1 ≤ j ≤ n − 1. The next result extends Theorem 3.1 of [12] to cover the case when a(x) satisfies (9).
If a, f satisfy (9) and (3), then (1) possesses a positive or negative solution if and only if λ Moreover, Proof. The proof follows mutatis mutandis from the proof of [12,Th. 3.1], except for the construction of the supersolution Φ µ in the proof of Proposition 3.1 of [12]. This construction must be modified slightly to cover the more general situation when (9) occurs. Under condition (9), given µ ∈Λ 1 , choose a sufficiently small δ > 0 such that Then, for each 1 ≤ j ≤ n, the function is positive in (α j − δ, β j + δ) and it satisfies Consider the function Adapting the argument of the proof of Proposition 3.1 of [12], it is easily seen that, for sufficiently large κ > 1, the function Φ µ = κφ δ is a supersolution of (1) for every λ ∈ (( π 2L ) 2 , µ]. Therefore, with this observation, the proof of [12,Th. 3.1] implies Theorem 2.1.
As in [12], the uniqueness of the positive and the negative solutions can be derived from the next extension of Theorem 2.2 of [12].
Suppose a, f satisfy (9) and (3) and λ > 0. Let u, v be solutions of (1) with u > 0 in (−L, L). Then, v ≤ u, i.e., u is the maximal solution of (1) and is the unique solution of (1) that is positive in (−L, L). Similarly, any negative solution of (1) is a minimal solution of (1) and is unique.
According to Theorem 2.2, for each λ ∈Λ 1 , the problem (1) has a unique positive solution, u λ , and a unique negative solution, v λ . Therefore, as a byproduct of Theorem 2.1, the next result holds.

Corollary 1.
Suppose a, f satisfy (9) and (3). Then, Similarly, the same statements hold if C + 1 and u λ are replaced by C − 1 and v λ . Actually, the following stronger statement which is a counterpart of Corollary 3.2 of [12] holds.
The next result we require is [12,Cor. 3.6] which was already stated in the form we need: Lastly, we will also need a version of Theorem 3.2 of [12].
Suppose a, f satisfy (9), (3), λ > 0 and Then, there are n + 1 functions, j (x), j ∈ {0, 1, ..., n}, such that Moreover, 0 is the minimal positive solution of the singular problem n is the minimal positive solution of and, for every 1 ≤ j ≤ n − 1, j is the minimal positive solution of Proof. By Corollary 2, u λ (x) has a limit, possibly ∞, as λ → ( π h ) 2 for each x ∈ [−L, L] and (11) holds for x ∈ ∪ n j=1 [α j , β j ]. To verify (16) for the remaining values of x, we adapt the argument of the proof of [12,Th. 3.2]. Consider the problems: and, for every 1 ≤ j ≤ n − 1, According to [8], or [9, Ch. 4], for every M > 0 and λ ∈ R, each of these problems possesses a unique positive solution, denoted by 0,M,λ , n,M,λ and j,M,λ , respectively. Moreover, for all 1 ≤ j ≤ n − 1, where 0,λ , n,λ and j,λ , 1 ≤ j ≤ n − 1, are the minimal positive solutions of (20), (21) and (22) with M = ∞, convergence being monotone in M for each x. The rest of the proof now readily follows from the argument in the proof of [12, Th. 3.2], so we omit any further technical details here. (15) holds if e.g. f (ξ) = |ξ| p for all ξ ∈ R and some p > 1. See e.g., [8], [9]). (ii) Let U be the limit function defined in (16). It is the minimal metasolution of (1). By a metasolution, of (1) 3. An example. In this section, an example will be given for which, under certain assumptions, the set of solutions with one node possesses an unbounded component, This does not exclude the possibility that either T 2 = C 2 , or Choose a weight function a = a 0 , symmetric about 0, such that, for some con- Actually, by Theorems 2.1 and 2.2, the problem has a unique positive solution for each λ ∈Λ 2 . So, due to the evenness of f , taking the unique positive solution of (24) and reflecting it as an odd function about x = 0 yields a one node solution of (1). In this fashion, for every λ ∈Λ 2 , (1) has a solution, (λ, u), with a node at 0 and and for ε ≥ 0, set a = a ε ≡ a 0 + εb.
We further assume a and b are such that a ε has a finite number of local maxima in (β 3 , α 4 ) for 0 ≤ ε ≤ ε 0 .
The next lemma provides us with the upper bound for z that will be required later.
Let V n,ε denote the unique negative solution of Since z(λ n , ε) < κ, the function U λn,ε provides us with a subsolution of (29) and consequently, by the maximum principle, Let W n denote the unique negative solution of (30). By (3) and the maximum principle again, V n,ε ≤ W n in [κ, L] for all n ≥ 1.
To obtain T 2 , we would like to invoke Corollary 4 and take T 2 = T . As was mentioned earlier, for any σ ∈ Λ 1 , the existence of an interval of values of ε for which the first hypothesis of Corollary 4 holds seems to be a reasonable assumption. The second hypothesis involving T is a rather strong condition. The numerical experiments carried out in Section 4 as well as the experience of first two authors in computing solution branches in a variety of reaction-diffusion equations and systems indicate that this hypothesis is valid in several settings. Thus we make the Conjecture: Suppose for some σ ∈ Λ 1 , u [σ,1, ] is a one node solution satisfying condition (M). Let T denote the component of the set of one node solutions containing u [σ, 1, ] . Then [σ, (π/h) 2 ) ⊂ PT and whenever λ ∈ (σ, (π/h) 2 ), there is a u such that (λ, u) ∈ T and u satisfies condition (M).

4.
A numerical example. Throughout this section, we take L = 1, f (ξ) = ξ 2 for all ξ ∈ R, and for sufficiently small ε > 0. Figure 2 shows a plot of this function for the choice ε = 0.1. For sufficiently small ε > 0, (1) can be regarded as a perturbation from the weight function a 0 obtained for the special choice ε = 0, whose graph has been already plotted in Figure 1. Thus, in this example, α = α 4 = 0.6 and β = β 4 = 0.8.
To discretize (1), a pseudo-spectral method based on [15] is used. To compute the subsequent global bifurcation diagrams we have adapted the path-following techniques of [1], [5], [7] and [10]. In the next subsections, the results of a series of numerical experiments for 3 significative choices of ε, namely 0.1, 0.0037, and 0.0036 will be discussed. Then, in the final section, some illuminating conclusions will be drawn from them.

4.1.
Numerical results for the choice ε = 0.1. The left plot of Figure 5 shows the computed global bifurcation diagram for π 2 < λ < 60. Using the notation of [12], it consists of a principal curve, C + 2 , bifurcating from u = 0 at λ = π 2 , together with an additional isola, I + , separated from the principal branch. The existence of this isola is not predicted by the theory of [12]. Only solutions (λ, u) with u (0) > 0 have been computed. Due to the choice of f , solutions of (1) occur in pairs (λ, ±u), but the solutions (λ, u) with u (0) < 0 will not be be plotted in the following bifurcation diagrams. The plot on the right magnifies the most significant piece of the left one. It shows that a kind of so-called imperfect bifurcation occurs as ε is perturbed from ε = 0. Regarding them as steady states of the parabolic counterpart of (1), the solutions on the upper half branch of the isola are stable, while those on the lower half of the branch are unstable with a one-dimensional unstable manifold.
The principal branch, C + 2 , has been computed until λ = 750, where the continuation step shortened to 10 −3 , making the computation of the one node solutions for further values of λ a very lengthy task. The solutions on the isola have been computed until λ = 190, where the continuation step shortened to 10 −3 with a high computational cost. The numerical evidence shows that for the special choice ε = 0.1, C + 2 = D + 2 and the solutions on the isola blow up in (α, β) = (0.6, 0.8) as λ ↑ 25π 2 . This is one of the alternatives, (ii), of Theorem 4.2 of [12], presented here in Theorem 1.1. Figure 6 gives a series of solution plots along the primary branch as λ varies. As in all subsequent solution plots, Figure 6 plots x as abscissa and u(x) as ordinate. It shows how the node of these solutions moves from 0.5 until 0.7 to avoid the blowup of the solution at λ = 25π 2 . According to Theorem 2.3 of [12], these are the unique solutions with one node of (1) for each λ ∈ [25π 2 , 100π 2 ). As λ increases and approximates 100π 2 , these solutions have to blow up in (0.6, 0.8) \ {0.7}, making their computation very costly due to the small size of the path-following step.  The left half of Figure 8 shows a series of solution plots along the isola for a further range of values of λ, 70 < λ < 140. These plots show how the nodes of the solutions on the lower branch, in blue, tend to z = 0, whereas the nodes of the solutions on the upper branch, in red, tend to z = −0.3. As a by-product, according to [12], all these solutions must blow up at λ = 25π 2 . The shape of the solutions on the isola for further values of λ, 140 < λ < 25π 2 , plotted in the right half of Figure  8 is very similar to those shown in the left half of the figure, although the bigger λ is, the more negative the solution becomes in (0.6, 0.8). As a consequence, if we wish to see their entire graphs in (α, β) = (0.6, 0.8), their plots should be rather flat in (−1, 0.6) and (0.8, 1). In Figure 9, the nodes of all the one node solutions of (1) that we have computed for λ ≤ 180 are plotted.  As predicted by [12], the node of the solutions on C + 2 = D + 2 tends to 0.7 as λ ↑ 100π 2 , while it tends to 0.0 along the upper half branch of the isola, and to −0.3 along the lower one. This explains why C + 2 can be continued until 100π 2 , the second eigenvalue of −d 2 /dx 2 on (0.6, 0.8) under 0 boundary conditions, whereas the isola blows up at λ = ( π h ) 2 .
As ε decreased from 0.1 down to 0.001, the previous isola persisted in approximating the principal branch. Figure 10 shows a magnification for the interval of λ's where they are close. But this isola is not important in our next numerical experiment. 19 19  Figure 10. A zoom of the bifurcation diagram for ε = 0.001

4.2.
Numerical results for the choice ε = 0.0037. As in the previous case our numerical computations strongly suggest that C + 2 = D + 2 also occurs in this case. Figure 11 shows the computed bifurcation diagram.  Figure 11. Two significant components of the bifurcation diagram Again C + 2 consists of linearly unstable one node solutions with one-dimensional unstable manifolds, except for a S-shaped piece around λ = 19.71, where the curve exhibited two consecutive turning points. This piece of C + 2 has been magnified in the left half of Figure 12. Indeed, a strong zoom in Figure 11 around λ = 19.71 reveals that there is a sort of protuberance in the bifurcation diagram of C + 2 around this value of the parameter. The left plot of Figure 12 shows a magnification of the strong zoom of Figure 11. The portion of the curve linking these turning points is filled in by linearly stable one node solutions. Throughout the rest of this paper, continuous lines are filled in by stable solutions and plotted in red while dashed lines are plotted in blue and are filled in by unstable solutions with one-dimensional unstable manifolds. 19   The second plot of Figure 12 shows a magnification of the unique turning point of the computed isola, J + , which cannot easily be detected in the global diagram plotted in Figure 11 because of the scale chosen to plot these two components. The relevant components in the next discussion are the ones plotted in Figures 11 and  12. Figure 13 uses a continuous line to plot of the zeroes of the one node solution along C + 2 = D + 2 and a doted-dashed line for the isola J + .
The left picture shows a series of solutions from λ = π 2 until λ = 60. The solutions for values λ ∼ π 2 are very small and have a shape similar to that of the second eigenfunction of −d 2 /dx 2 in (0, 1) under 0 Dirichlet conditions. In particular, their zeroes are very close to 0.0. As λ increases from π 2 , their zeroes move toward the right approximating 0.5 as λ tends to 60, where the magnification shown in the right half of Figure 14 was taken. In the second half of Figure 15, the graphs of a series of one node solutions (λ, u) with u (0) > 0 on C + 2 has been plotted for a range of values of λ between π 2 and 700. Due to the the scale used, the plots for the smaller values of λ cannot be distinguished. According to Theorem 2.3 of [12], for every λ ∈ [25π 2 , 100π 2 ), (1) possesses a unique one node solution u λ := u [λ, 1,ε] with u λ (0) > 0. Moreover, by Remark 5.7 of [12], where stands for the unique positive solution of the singular problem x ∈ (0, 0.6), u(0) = 0, u(0.6) = +∞, and m is the unique negative solution of the singular problem Note that the right side plot of Figure 15 illustrates this limiting behavior. Figure 16 shows  Figure 17 shows a series of solution plots around λ = 19.71 in order to see what's happens to the solutions in the region between the two turning points along the component C + 2 . Essentially, these turning points facilitate the formation of the third hump of the solutions, speeding the movement of the node towards the right so that it can reach 0.6 before λ reaches 25π 2 . The right picture focuses on a piece of the plots of all solutions of C + 2 = D + 2 allowing the reader to see the convergence of their interior nodes to 0.7 as λ ↑ 100π 2 . Indeed, the nodes approximate the right end of the x-axis, at 0.7.

4.3.
Numerical results for the choice ε = 0.0036. As ε decreases from ε = 0.0037 down to ε = 0.0036, the previous situation changes drastically. Indeed, the bifurcation diagram of Figure 11 perturbs into the bifurcation diagram plotted in Figure 19.
Although at first glance the bifurcation diagrams of Figures 11 and 19 look very similar, they are rather different. Indeed, as is illustrated in Figure 20, the bifurcation diagram consists of two different components, C + 2 and D + 2 . Now, the component C + 2 blows up at λ = 25π 2 and D + 2 blows up at 25π 2 and 100π 2 , in complete agreement with Theorem 1.1(ii).   Figure 19 In this case the zeroes of the one node solutions on C + 2 do not reach the interval (α, β) = (0.6, 0.8) before λ attains the critical value of 25π 2 . This explains why C + 2 blows up at λ = 25π 2 for this value of ε. Figure 21 shows the plot of the zeroes of the one node solutions of (1) along the components C + 2 using a solid line, and D + 2 with a dashed-dotted line. It should be compared with the diagram plotted in Figure 13. Figure 22 shows two consecutive magnifications around λ = 64 of the curves of zeroes plotted in Figure 21. They should be compared with the corresponding magnifications in Figure 14. The nodes of the solutions along C + 2 grow and approach 0.5 but just before reaching this value, around 0.495, they begin to decrease for all further values of λ. Consequently, the solutions on C + 2 must blow up in (α, β) = (0.6, 0.8) as λ ↑ 25π 2 . The limiting behavior of these solutions as λ ↑ 25π 2 can be described by the same large solutions as in case ε = 0.0037 along the isola J + . The plots of the solutions on the components C + 2 for π 2 < λ < 60 are very similar to the plots of the corresponding solutions for ε = 0.0037. Similarly, the plots of the solutions on C + 2 for 60 < λ < 25π 2 are very similar to the plots of the solutions along the lower half branch of the isola J + computed for ε = 0.0037. Similarly, the plots of the solutions on the lower half branch of D + 2 for 65 < λ < 100π 2 are perturbations of the ones computed along C + 2 for this range of values of λ for ε = 0.0037. Naturally, the plots of the solutions on the upper half branch of D + 2 , those in the continuous line in Figure 19, are perturbations of those of the stable one node solutions on the upper half branch of the isola J + 2 for ε = 0.0037. Consequently, we will not carry out any further solution plots here.
5. Discussion and conclusions. All our numerical experiments showed that, for every ε ≥ 0.0037, C + 2 = D + 2 , whereas C + 2 ∩ D + 2 = ∅ if ε ≤ 0.0036. Consequently, the numerical experiments strongly suggest that the situation described in Theorem 1.1 (i) occurs if ε ≥ 0.0037, whereas the situation described by Theorem 1.1 (ii) occurs under condition ε ≤ 0.0036. One suspects that by continuous dependence, there should exist a ε * > 0 satisfying 0.0036 < ε * < 0.0037 where C + 2 = D + 2 blows up, simultaneously, at 25π 2 = ( π β−a ) 2 and 100π 2 = ( 2π β−α ) 2 , which also fits into the abstract setting of Theorem 1.1 (i). The bifurcation diagram for ε = ε * should be very similar to the one plotted in Figures 11 and 19, but with a unique component. Likewise, the zeroes diagram should be similar to the one shown in the left plots of Figures 14 and 22. All these conjectures, which are supported by the numerical experiments we have carried out, are based on the fact that the equation −u = λu − a ε (x)u 3 x ∈ (0.5, 1), u(0.5) = u(1) = 0, admits a unique negative solution. Consequently, if the node approximates 0.5 as ε → ε * , the solutions should converge to a common profile which should be a branching point to four different solution branches. Naturally, as ε perturbs from ε * , the component spreads out into two components, whose structure depends on whether ε > ε * , or ε < ε * .
Thus this example indicates how all the cases enumerated in our Theorem 1.1, which is Theorem 4.2 of [12], can indeed occur by choosing an appropriate ε. Moreover, both the heuristics and the computations indicate these models possess at least an additional unbounded isola of solutions, I + , with the turning point around 19, which blows up to infinity at λ = 25π 2 . Figure 1 shows it for ε = 0.1 and Figure  10 shows a piece of it around 19 for ε = 0.001.
Since the previous examples possesses at least three solutions for each λ ∈ (65, 25π 2 ), the uniqueness theorem established by Theorem 2.3 of [12] is optimal in that for λ / ∈ (25π 2 , 100π 2 ), the uniqueness of the one node solution cannot be guaranteed.
The occurrence of isolas of positive solutions for a class of stepwise constant weight functions, a(x), that change sign has been demonstrated recently in [13]. The emergence of isolas in the context of the current paper, despite that a ≥ 0 in (1), is attributable to the fact that, as far as nodal solutions are concerned, (1) behaves like an indefinite problem in the sense that f (u(x)) changes sign if u(x) changes sign. This explains why it is a considerable challenge to ascertain whether or not C + 2 = D + 2 . It should be noted that in our spatially heterogeneous context, we cannot invoke the nice phase portrait techniques of [13].