Electronic Journal of Differential Equations, Vol. 2025 (2025), No. 92, pp. 1–19.

ISSN: 1072-6691. URL: https://ejde.math.txstate.edu, https://ejde.math.unt.edu

DOI: 10.58997/ejde.2025.92

VARIATIONAL AND NUMERICAL ASPECTS OF A SYSTEM OF ODES

WITH CONCAVE-CONVEX NONLINEARITIES

OSCAR AGUDELO, GABRIELA HOLUBOVÁ, MARTIN KUDLÁČ

Abstract. We study a Hamiltonian system of ordinary differential equations under Dirichlet

boundary conditions. The system features a mixed (concave-convex) power nonlinearity with a

positive parameter λ. We show multiplicity of nonnegative solutions for a range of the parameter
λ and discuss the regularity and symmetry of nonnegative solutions. Besides, we present a

numerical strategy aiming at the exploration of the optimal range of λ for which multiplicity

of positive solutions holds. The numerical experiments are based on the Poincaré-Miranda
Theorem and the shooting method, which have been lesser explored in the context of multiple

positive solutions of systems of ODEs. Our work has been motivated by the results in Ambrosetti

et al. in [4] and Moreira dos Santos in [18].

1. Introduction

We study the system of ordinary differential equations (ODEs for short)

−u′′ = λ|v|r−1v + |v|p−1v in (0, 1),

−v′′ = |u|q−1u in (0, 1)
(1.1)

with Dirichlet boundary conditions

u(0) = u(1) = 0 and v(0) = v(1) = 0. (1.2)

Here λ is a positive parameter and the exponents p, q and r are assumed to satisfy

0 < r <
1

q
and p > max

{
1,

1

q

}
. (1.3)

Systems of equations such as those in (1.1) appear in the study of population dynamics, fluid
dynamics and stellar structure in astrophysics. Related systems of equations have been explored
since the nineteenth century to study for instance the density of a gas sphere (see e.g. [13, 24]).

System (1.1) features concave-convex polynomial nonlinearities. When 0 < r < 1 < p, the
nonlinearity in the first equation of (1.1) is concave near the origin and convex at infinity. When
0 < q < 1 < r, the nonlinearity in the first equation of (1.1) stays convex while concavity appears
in the second equation of (1.1). We refer the reader to [11] for the study of qualitative properties
and the classification of radial solutions to some problems related to (1.1), which appear in the
modeling of thermal structures of static configurations, and in which the transfer of energy takes
place entirely by thermal conduction. The nonlinearities in [11] involve concave and convex terms,
accounting for heat radiation and heat generation, respectively.

In [4], the authors study existence, nonexistence and multiplicity of nonnegative solutions of
the single equation

−v′′ = λ |v|r−1v + |v|p−1 v in (0, 1) (1.4)

2020 Mathematics Subject Classification. 34A34, 34B08, 34B18, 35J35.
Key words and phrases. Hamiltonian system of odes; concave and convex nonlinearities; minimization theorem;

mountain pass theorem; shooting method; moving planes method.
©2025. This work is licensed under a CC BY 4.0 license.
Submitted January 12, 2025. Published October 3, 2025.

1



2 O. AGUDELO, G. HOLUBOVÁ, M. KUDLÁČ EJDE-2025/92

with Dirichlet boundary conditions
v(0) = v(1) = 0 (1.5)

and with respect to the nonnegative parameter λ, assuming that 0 < r < 1 < p < +∞. Actually,
the study in [4] treats mainly the higher dimensional case, where technical issues related to the
non-compactness of certain Sobolev embeddings arise. In [18], the author studies system (1.1)-
(1.2) exclusively in higher dimensions. The spirit of the results in [18] is similar to the one in
[4]. In [18], the range of values λ for which existence and multiplicity of positive solutions of
(1.1)-(1.2) hold, is not optimal and we remark that the techniques in [18] differ from the ones in
[4]. In [2], the authors discuss the existence of minimal solutions for a related Hamiltonian system
of equations. Nonexistence of solutions is also discussed. We refer the reader to the surveys [16]
and [28] and references therein for a description of some of the known results related to systems
of differential equations and the techniques used to treat them.

Motivated mainly by [4] and [18], we explore in dimension one the existence and multiplicity
of nonnegative classical solutions of (1.1)-(1.2) with respect to the parameter λ. By a classical
nonnegative solution of (1.1)-(1.2) we mean a pair (u, v) of functions such that u, v ∈ C2(0, 1) ∩
C[0, 1], u, v ≥ 0 in (0, 1) and (1.1) and (1.2) being satisfied point-wise. Our main result reads as
follows.

Theorem 1.1. There exists a positive constant λ0 such that for any λ ∈ (0, λ0), system (1.1)-(1.2)
possesses at least two distinct nontrivial nonnegative classical solutions.

We now explain briefly the strategy of the proof of Theorem 1.1. Based on the method of
reduction by inversion (see e.g. [14, 23]), and proceeding as in [18], we reformulate (1.1)-(1.2) as
a single fourth order ODE with Navier boundary conditions (see BVP (3.3)-(3.4)). This fourth
order BVP has variational structure (for the corresponding energy functional see (3.6)). Following
partly the approach in [4], we apply the Ambrosetti-Rabinowitz Mountain Pass Theorem and a
local minimization argument to the associated energy functional to prove the existence of two
distinct solutions for sufficiently small λ (see Propositions 4.3 and 5.2, respectively).

We remark that in contrast with the fibering method used in [18], this variational approach
easily adapts to study existence of nonnegative solutions of BVP of the type (1.1)-(1.2) with more
general nonlinearities. For instance, one may consider nonlinearities with similar mixed polynomial
growth, but such that this feature does not come from autonomous terms. For the sake of clarity
and emphasis of the main ideas in our exposition, we work directly with the BVP (1.1)-(1.2).

In addition to discussing existence and multiplicity of classical solutions, we prove the regularity
of weak solutions of the associated fourth order BVP (see BVP (3.3)-(3.4) and Proposition 3.1) and
provide a lower estimate of the optimal value of λ0. We finish our theoretical discussion by using
the method of moving planes to prove symmetry of the solutions of (1.1)-(1.2) (see Proposition
8.1).

Although we use quite standard functional analytical tools, the nonlinear terms in (1.1) require
a fine and careful treatment. Due to the techniques and restrictions of the Sobolev embeddings
used in [18], the results concerned with existence of solutions for the higher-dimensional version
of system (1.1) do not directly translate into the one-dimensional setting. Nonetheless, as in [18],
the variational techniques used in this work do not provide the optimal quantitative information
concerning the range of values λ for which existence and multiplicity of nonnegative solutions hold
true. Its analytic description (depending on the parameters p, q, r) remains to be open. In this
regard, recently in [3] the authors explore the theoretical aspects of the questions posed herein,
but for a slightly more general system of equations.

Thus, we further explore this matter using a numerical implementation motivated by a com-
bination of the standard shooting method and the heuristics of the Poincaré-Miranda Theorem.
One of the first numerical explorations in this direction is carried out in [27], where a basic im-
plementation of the shooting method is developed to study BVPs mainly for linear systems of
equations. On the other hand, relying on different versions of the Intermediate Value Theorem,
the shooting method has been greatly exploited to study multiplicity of sign-changing solutions of
BVPs in different contexts. In this regard, we refer the reader to [8, 9, 10, 12, 15] and references
therein.



EJDE-2025/92 ODES WITH CONCAVE-CONVEX NONLINEARITIES 3

To the best of our knowledge, this strategy has been rarely explored to study multiplicity of
positive solutions of systems of the type (1.1)-(1.2) and hence we believe the numerical illustration
presented in Section 9 is a novel, insightful and versatile approach that easily adapts to more
general systems of ODEs.

In particular, in line with the statement of Theorem 1.1, we obtain a bifurcation diagram
showcasing the dependence of the L∞(0, 1)-norm of the v-component of a solution of (1.1)-(1.2)
with respect to the parameter λ (see Figure 1).

λ

m
a
x
x
∈
[0
,1
]
v
(x
)

Figure 1. Bifurcation diagram for the parameter λ. Here, p = 3, q = 3/2,
r = 1/3

Borrowing some terminology from the Bifurcation Theory, our results suggest the following:
there exists λbif > 0 such that for λ ∈ [0, λbif) there are at least two solutions – a stable one vλ
and the unstable one vλ. At the point of bifurcation λ = λbif , the solutions vλ and vλ coincide.
We remark that, as shown in [2, Theorem 1.3], when 0 < r < 1

q < 1 < p, there exists Λ0 ∈ (0,∞)

such that for any λ > Λ0 system (1.1) has no nontrivial nonnegative solutions. It is expected that
the smallest of such Λ0 coincides with λbif .

This article is organized as follows. Section 2 introduces preliminary results and notion required
in subsequent sections. In Section 3 we set the functional analytic framework to study problem
(1.1)-(1.2). This section also discusses the regularity of solutions. To carry out the proof of the
main result, Sections 4 and 5 provide existence of two weak solutions of the fourth order BVP (see
(3.3)-(3.4)) via the Mountain Pass Theorem and a local minimization argument for the associated
energy functional (see (3.6)), respectively. Section 6 deals with compactness properties of approx-
imating sequences of solutions, namely, the Palais-Smale condition. The proof of Theorem 1.1 is
carried out in Section 7. Section 8 is devoted to the proof of Proposition 8.1 describing symmetry
of the solutions. The last section presents various numerical illustrations with detailed description
of our numerical strategy.

2. Preliminaries

Let s ∈ [1,∞] and let I be a bounded open interval. In what follows, Ls(I) denotes the Lebesgue
space of measurable functions w : I → R endowed with the norm

∥w∥Ls(I) :=

{( ∫
I
|w|sds

)1/s
, s ∈ [1,+∞),

ess supx∈I |w(x)|, s = +∞.

When I = (0, 1) and for k ∈ N, we also consider the Sobolev space W k,s(0, 1), which consists
of all functions u ∈ Ls(0, 1) such that for any i ∈ {1, . . . , k} the i−th weak derivative of u, u(i),



4 O. AGUDELO, G. HOLUBOVÁ, M. KUDLÁČ EJDE-2025/92

belongs to Ls(0, 1). The space W k,s(0, 1) is endowed with the norm

∥u∥Wk,s(0,1) := ∥u∥Ls(0,1) +

k∑
i=1

∥u(i)∥Ls(0,1).

Recall that the Sobolev space W k,s(0, 1) is a Banach space. Furthermore, W k,s(0, 1) is separable
for s ∈ [1,∞) and reflexive for s ∈ (1,∞) (see [1, Theorems 3.3 and 3.6, p. 60–61]).

Let s > 1 and let W 1,s
0 (0, 1) denote the closure of C1

c (0, 1) in W 1,s(0, 1) with respect to the

norm ∥.∥W 1,s(0,1). A convenient description of the space W 1,s
0 (0, 1) is the following (see [7, Th.

8.12, p. 217]): for any u ∈ W 1,s(0, 1),

u ∈ W 1,s
0 (0, 1) if and only if u = 0 on ∂I.

Lemma 2.1 (Morrey’s inequality revisited). If u ∈ W 2,s(0, 1) ∩W 1,s
0 (0, 1), then

∥u∥L∞(0,1) ≤
1

2
∥u′′∥Ls(0,1).

Proof. First, notice that if w ∈ L∞(0, 1), then for any x, y ∈ (0, 1) with x < y, we estimate

∥w∥Ls(x,y) =
(∫ y

x

|w|sds
)1/s

≤ (y − x)1/s∥w∥L∞(x,y). (2.1)

Also, if w ∈ W 1,s(0, 1), then for any x, y ∈ [0, 1] with x < y, w ∈ W 1,s(x, y). Even more, from the
Morrey’s inequality ([20, Th. 4, p. 280]),

|w(y)− w(x)| ≤ (y − x)1−
1
s ∥w′∥Ls(x,y). (2.2)

Finally, if w ∈ W 1,σ
0 (0, 1) with σ ≥ 2, then, for any x ∈ (0, 1), (2.2) implies

|w(x)|σ ≤ xσ−1∥w′∥σLσ(0,x) and |w(x)|σ ≤ (1− x)σ−1∥w′∥σLσ(x,1).

Multiplying the left inequality by (1 − x)σ−1 and the right inequality by xσ−1 and adding these
terms yield (

(1− x)σ−1 + xσ−1
)
|w(x)|σ ≤ xσ−1(1− x)σ−1∥w′∥σLσ(0,1).

Applying Jensen and Arithmetic-Geometric mean inequalities, we obtain

|w(x)| ≤ 1

2
∥w′∥Lσ(0,1). (2.3)

Now, let u ∈ W 2,s(0, 1)∩W 1,s
0 (0, 1), s > 1. Then u ∈ W 1,σ

0 (0, 1) with any σ ≥ 2 and thus (2.2),
(2.3), and Rolle’s Theorem (see [30, p. 215]) yield

∥u∥L∞(0,1) ≤
1

2
∥u′∥Lσ(0,1) ≤

1

2
∥u′∥L∞(0,1) ≤

1

2
∥u′′∥Ls(0,1),

in other words the statement of this lemma. □

Next, let us fix q > 0 and introduce

X := W 2, q+1
q (0, 1) ∩W

1, q+1
q

0 (0, 1). (2.4)

Observe that u ∈ X implies that u ∈ C1, 1
q+1 [0, 1], u′′ ∈ L

q+1
q (0, 1) and u(0) = u(1) = 0.

It is readily verified using Lemma 2.1 that

∥v∥X :=
(∫ 1

0

|v′′(x)|
q+1
q dx

) 1
q+1

is a norm in X. From now on (unless stated otherwise), we endow X with this norm.

Remark 2.2. Notice that Lemma 2.1 implies the following. If u ∈ X, then for any s ∈ (1,∞],
u ∈ Ls(0, 1) and

∥u∥Ls(0,1) ≤ ∥u∥L∞(0,1) ≤
1

2
∥u∥X . (2.5)

Moreover, the linear embedding

i : X → Ls(0, 1), i(u) := u (2.6)



EJDE-2025/92 ODES WITH CONCAVE-CONVEX NONLINEARITIES 5

is continuous in X and ∥i∥ ≤ 1
2 .

Lemma 2.3. The normed space X has the following properties:

(a) the norms ∥ · ∥X and ∥ · ∥
W

2,
q+1
q (0,1)

are equivalent in X,

(b) the space X is Banach,
(c) the space X is reflexive,
(d) the space X is compactly embedded into C1[0, 1].

Claims (a)–(c) are thoroughly proved in [25], the last property of X is obtained by iterating of

compact embedding of W 1, q+1
q (0, 1) into C[0, 1] (see [7, Th. 8.2, p. 204]).

3. Functional analytic setting

In this section we discuss several formulations of (1.1)-(1.2) and concepts of its solutions. Recall
that p, q and r satisfy (1.3).

First of all, to capture nonnegative solutions, instead of working directly with (1.1)-(1.2), we
consider the system

−u′′ = λ(v+)
r + (v+)

p in (0, 1),

−v′′ = |u|q−1u in (0, 1),

u(0) = u(1) = v(0) = v(1) = 0,

(3.1)

where t+ := max{t, 0} for t ∈ R. Indeed, let u, v ∈ C2(0, 1)∩C[0, 1] be such that the pair (u, v) is
a classical solution of (3.1), that is, all the identities in (3.1) hold in point-wise sense. First notice
that if the pair (u, v) is not the trivial vector function (0, 0), then both u and v are nontrivial.
Next, from the first equation in (3.1), u is concave in [0, 1]. Moreover, since u(0) = u(1) = 0, u
is nonnegative in [0, 1]. Arguing in a similar manner with the second equation, using that u ≥ 0
in (0, 1), the same holds for v. In conclusion, any classical nontrivial solution (u, v) of (3.1) is
such that u and v are positive in (0, 1) and hence it is also a solution of (1.1)-(1.2). Clearly, any
classical solution (u, v) of (1.1) with u ≥ 0 and v ≥ 0 in [0, 1] solves also (3.1).

Also, from the second equation in (3.1) we have

u = −|v′′|
1
q−1v′′ ∈ C2(0, 1). (3.2)

Plugging (3.2) into the first equation in (3.1), system (1.1)-(1.2) reduces to the fourth order
equation

d2

dx2

(
|v′′|

1
q−1v′′

)
= λvr+ + vp+ in (0, 1) (3.3)

with the Navier boundary conditions

v(0) = v(1) = 0 and v′′(0) = v′′(1) = 0. (3.4)

A classical solution of (3.3)-(3.4) is a function v ∈ C2[0, 1] such that |v′′|
1
q−1v′′ ∈ C2(0, 1)∩C[0, 1]

and (3.3)-(3.4) are satisfied point-wise. Observe that nonnegative classical solutions of system
(1.1)-(1.2) are in correspondence with classical solutions of BVP (3.3)-(3.4).

Besides the classical setting for solutions, we will use the concept of weak solutions. We consider
the space X defined in (2.4). By a weak solution to (3.3)-(3.4) we mean a function v ∈ X such
that for any φ ∈ X, ∫ 1

0

|v′′|
1
q−1v′′φ′′ dx =

∫ 1

0

(
λvr+ + vp+

)
φdx. (3.5)

In the following statements we show that weak and classical solutions of BVP (3.3)-(3.4) coincide.

Proposition 3.1. If v ∈ X is a weak solution of (3.3)-(3.4), then v is also a classical solution of
(3.3)-(3.4).

Proof. Let v ∈ X be a weak solution of (3.3)–(3.4). We write w := |v′′|
1
q−1v′′ and h := λvr+ + vp+

a.e. in (0, 1). Observe that w ∈ Lq+1(0, 1), h ∈ C[0, 1] and from (3.5) for any φ ∈ X∫ 1

0

wφ′′ dx =

∫ 1

0

hφdx.



6 O. AGUDELO, G. HOLUBOVÁ, M. KUDLÁČ EJDE-2025/92

Let ω ∈ C2[0, 1] solve the BVP

ω′′ = h in (0, 1), ω(0) = ω(1) = 0.

We prove that w = ω a.e. in (0, 1). Once this is proven, we would conclude that |v′′|
1
q−1v′′ ∈

C2[0, 1].
To prove the claim, we write ϑ := |w − ω|q−1(w − ω) a.e. in (0, 1). Since w ∈ Lq+1(0, 1), we

have ϑ ∈ L
q+1
q (0, 1). Let φ ∈ W 2, q+1

q (0, 1) solve the BVP

φ′′ = ϑ in (0, 1), φ(0) = φ(1) = 0.

Observe that φ ∈ C1[0, 1] and from the boundary conditions for φ we find that φ ∈ X. Using φ
as a test function in (3.5), we obtain∫ 1

0

|w − ω|q+1dx =

∫ 1

0

(w − ω)ϑ dx (definition of ϑ)

=

∫ 1

0

(w − ω)φ′′dx (choice of φ)

=

∫ 1

0

wφ′′dx−
∫ 1

0

ω′′φdx (integration by parts)

=

∫ 1

0

hφdx−
∫ 1

0

hφdx = 0.

Thus, w = ω a.e. in (0, 1) and this proves the claim. Since v′′ = |w|q−1w in [0, 1], then v′′ ∈ C[0, 1].
Consequently, v ∈ C2[0, 1]. This completes the proof. □

For an alternative proof of Proposition 3.1 (under slightly stronger assumptions on the exponent
q, motivated by the arguments presented in [21]), we refer the reader to [25]. The following
corollary is a direct consequence of Proposition 3.1.

Corollary 3.2. Let v ∈ X \ {0} be a nonnegative (weak) solution of (3.3)-(3.4) and set u =

−|v′′|
1
q−1v′′ a.e. in (0, 1). Then u, v ∈ C2[0, 1], u, v > 0 in (0, 1) and the pair (u, v) is a classical

solution of (1.1)-(1.2).

Summarizing the above results, classical nonnegative nontrivial solutions to (1.1)-(1.2) are in
correspondence with the nontrivial weak solutions of (3.3)-(3.4). As we have already announced,
we will use a variational approach to find them.

Let us consider the energy functional J : X → R defined by

J(v) :=
q

q + 1

∫ 1

0

|v′′|
q+1
q dx− 1

r + 1

∫ 1

0

λvr+1
+ dx− 1

p+ 1

∫ 1

0

vp+1
+ dx. (3.6)

From Remark 2.2, J is well-defined. It is also standard to verify that J ∈ C1(X) with

DJ(v)φ =

∫ 1

0

|v′′|
1
q−1v′′φ′′dx−

∫ 1

0

(
λvr+ + vp+

)
φdx for v, φ ∈ X. (3.7)

Thus, DJ(v) = 0 if and only if v is a weak solution of (3.3)–(3.4). We refer the reader to [25] for
details.

For convenience, we write J(v) as

J(v) =
q

q + 1
∥v∥

q+1
q

X − λ

r + 1
∥v+∥r+1

Lr+1(0,1) −
1

p+ 1
∥v+∥p+1

Lp+1(0,1). (3.8)

In the following sections, we examine the critical points of J .



EJDE-2025/92 ODES WITH CONCAVE-CONVEX NONLINEARITIES 7

4. Mountain pass solution

In this section, we use the Mountain Pass Theorem to show that for λ > 0 small, (3.3)–(3.4)
has at least one solution. In what follows, X∗ denotes the topological dual space of X with the
topology induced by the norm in X. Recall that p, q, r satisfy (1.3).

First, we remark that the functional J satisfies the Palais-Smale condition (PS-condition for
short) in X, i.e., for any c ∈ R and for any sequence {un}n ⊂ X such that

J(un) → c and ∥DJ(un)∥X∗ → 0, (4.1)

there exists a subsequence {unk
}k ⊂ {un}n that converges strongly in X (see e.g. [19]). The

detailed proof of this fact is postponed until Section 6.
Next, we verify that the functional J has the Mountain Pass geometry. That is the scope of

the following lemmas.

Lemma 4.1. There exist positive constants T and λ0 = λ0(T, p, q, r) such that for all λ ∈ (0, λ0)
there exists C = C(T, p, q, r, λ) so that for any v ∈ X satisfying ∥v∥X = T , J(v) ≥ C > 0.

Proof. Observe that for any v ∈ X, v+ ≤ |v| in (0, 1). Using (3.8) and (2.5) from Remark 2.2 we
estimate

J(v) ≥ q

q + 1
∥v∥

q+1
q

X − λ

2r+1(r + 1)
∥v∥r+1

X − 1

2p+1(p+ 1)
∥v∥p+1

X

≥ ∥v∥r+1
X

( q

q + 1
∥v∥

1
q−r

X − λ

2r+1(r + 1)
− 1

2p+1(p+ 1)
∥v∥p−r

X

)
.

(4.2)

Consider the function

h(t) :=
q

q + 1
t
1
q−r − 1

2p+1(p+ 1)
tp−r, t ≥ 0.

A direct calculation yields that

h′(t) = t
1
q−r−1

(1− qr

q + 1
− p− r

2p+1(p+ 1)
tp−

1
q

)
for t > 0.

Using (1.3) we find that the only two possible points of local extrema of h(t) are

t1 = 0 and t2 =
(2p+1(1− qr)(p+ 1)

(q + 1)(p− r)

) q
pq−1

.

Again, inequalities in (1.3) yield h(0) = 0, limt→∞ h(t) = −∞, and h′ is positive for t ≪ 1. Thus,
t2 is a point of global maximum of h and h(t2) > 0. We denote T := t2.

Let v ∈ X satisfy ∥v∥X = T . Inequality (4.2) for v now reads as

J(v) ≥ T r+1
(
h(T )− λ

2r+1(r + 1)

)
. (4.3)

We denote

λ0 := 2r+1(r + 1)
( q

q + 1
T

1
q−r − 1

2p+1(p+ 1)
T p−r

)
= 2r+1(r + 1)h(T ) (4.4)

and let λ ∈ (0, λ0). Setting C :=
(
1− λ

λ0

)
T r+1h(T ), we observe that

C = T r+1

(
q

q + 1
T

1
q−r − 1

2p+1(p+ 1)
T p−r − λ

2r+1(r + 1)

)
.

It follows from (4.3) and (4.4) that for v ∈ X with ∥v∥X = T ,

J(v) ≥ C > 0,

which proves the claim. □

In what follows, T > 0 and λ0 > 0 are the constants stated in Lemma 4.1. Let BT ⊂ X denote
the closed ball centered at the origin with radius T , that is,

BT = {u ∈ X : ∥u∥X ≤ T}.

Lemma 4.2. For any given λ > 0, there exists e ∈ X \BT such that J(e) < 0.



8 O. AGUDELO, G. HOLUBOVÁ, M. KUDLÁČ EJDE-2025/92

Proof. Let φ ∈ X be an arbitrary, but fixed function such that φ ∈ ∂BT and φ > 0 in (0, 1).
Using (3.8), for t positive

J(tφ) =
q

q + 1
t
q+1
q ∥φ∥

q+1
q

X − λ

r + 1
tr+1∥φ∥r+1

Lr+1(0,1) −
1

p+ 1
tp+1∥φ∥p+1

Lp+1(0,1), (4.5)

and from (1.3) we have

lim
t→∞

J(tφ) = −∞.

Therefore for a sufficiently large t > 1, the function e := tφ is such that e ∈ X \ BT and
J(e) < 0. □

Now we are in a position to prove the existence of a first solution to (3.3)–(3.4).

Proposition 4.3. Let λ0 > 0 be as in Lemma 4.1. For λ ∈ (0, λ0), there exists a nontrivial weak
solution v1 ∈ X of (3.3)–(3.4). Moreover, J(v1) > 0.

Proof. Observe that J ∈ C1(X,R) and satisfies the PS-condition (see Section 6). Also, J(0) = 0
and from Lemma 4.1, for λ ∈ (0, λ0), infv∈∂BT

J(v) ≥ C > 0. Finally, from Lemma 4.2, there
exists e ∈ X such that ∥e∥X > T and J(e) < 0. Let

c := inf
γ∈Γ

max
0≤s≤1

J(γ(s)),

where Γ := {γ ∈ C([0, 1];X) : γ(0) = 0 and γ(1) = e}. The Mountain Pass Theorem (see [19,
Theorem 6.4.5] with F = J and R = T ) yields the existence of v1 ∈ X such that

J(v1) = c and DJ(v1) = 0 in X∗. (4.6)

The definition of c yields that J(v1) ≥ C > 0. From (4.6), v1 is a nontrivial weak solution of
(3.3)–(3.4). This completes the proof. □

5. Local minimum solution

We recall again that p, q, r satisfy (1.3). In this part we show that for λ > 0 small enough,
(3.3)–(3.4) has a nontrivial second solution which is a local minimum for J near the trivial solution.

We use the notations and conventions introduced in Section 4. Let T , λ0 and C be as in the
statement of Lemma 4.1 and recall that BT is the closed ball in X centered at zero with radius T .

Lemma 5.1. Let λ ∈ (0, λ0). There exists ṽ ∈ X \ {0} such that ṽ ∈ intBT and J(ṽ) < 0.

Proof. Let φ ∈ X be an arbitrary, but fixed function such that φ ∈ ∂BT and φ > 0 in (0, 1).
Then, using (4.5), for t ∈ (0, 1) and λ > 0,

J(tφ) = tr+1
( q

q + 1
t
1
q−r∥φ∥

q+1
q

X − λ

r + 1
∥φ∥r+1

Lr+1(0,1) −
1

p+ 1
tp−r∥φ∥p+1

Lp+1(0,1)

)
.

Since 1
q > r > 0, λ ∈ (0, λ0) and since T > 0 does not depend on λ, we may select t = O(λ

q
1−qr )

small enough such that t ∈ (0, T ) and J(tφ) < 0. Finally, if we denote ṽ := t φ, then ṽ ∈ BT and
J(ṽ) < 0. This completes the proof. □

The existence of another solution is stated next.

Proposition 5.2. Let T > 0 and λ0 be as in Lemma 4.1. For all λ ∈ (0, λ0), there exists a
nontrivial weak solution v2 ∈ intBT of (3.3)–(3.4). In addition, J(v2) < 0.

Proof. Consider the minimization problem

mλ := inf
v∈BT

J(v). (5.1)

We claim that mλ is attained. To prove this we proceed as follows. Recall that X is reflexive (see
Lemma 2.3) and thus, using Kakutani’s Theorem (see [7, Theorem 3.17, p. 67]), BT is sequentially
weakly compact. It suffices to prove that J is sequentially weakly lower semicontinuous in BT .
Once this is proven, the existence of v2 ∈ BT such that J(v2) = mλ will follow from [29, Theorem
1.1], and the corresponding comments. In particular, mλ ∈ R. Also, Lemma 4.1 and Lemma



EJDE-2025/92 ODES WITH CONCAVE-CONVEX NONLINEARITIES 9

5.1 guarantee that v2 ∈ intBT with J(v2) < 0. Consequently, v2 is a nontrivial weak solution of
(3.3)–(3.4).

Now we prove the sequential weak lower semicontinuity of J in BT . Let {vn}n ⊂ BT be
arbitrary and such that vn ⇀ v weakly in X. Since BT is sequentially weakly compact, v ∈ BT .
We prove that

J(v) ≤ lim inf
n→+∞

J(vn).

Notice that X is continuously embedded in W 1, q+1
q (0, 1) and hence compactly embedded into

C[0, 1]. Thus, the sequence {vn}n converges uniformly to v in [0, 1]. In particular, (vn)+ → v+
uniformly in [0, 1]. Since the norm ∥ · ∥X is sequentially weakly lower semicontinuous in BT [7,
Prop. 3.5 (iii), p. 58], we estimate

lim inf
n→∞

J(vn) ≥
q

q + 1
∥v∥

q+1
q

X − λ

r + 1
lim

n→∞

(∫ 1

0

(vn)
r+1
+ dx

)
− 1

p+ 1
lim
n→∞

(∫ 1

0

(vn)
p+1
+ dx

)
=

q

q + 1
∥v∥

q+1
q

X − λ

r + 1

(∫ 1

0

vr+1
+ dx

)
− 1

p+ 1

(∫ 1

0

vp+1
+ dx

)
= J(v).

(5.2)

This proves the last claim and completes the proof. □

6. Palais-Smale condition

Lemma 6.1. The energy functional J satisfies PS-condition.

Proof. Let us consider a sequence {vn}n ⊂ X satisfying (4.1) with c ∈ R. This implies that there
exists K > 0 such that

J(vn) ≤ K ∀n ∈ N. (6.1)

It also implies that for any ε > 0 there exists n0 ∈ N such that for all n ∈ N with n ≥ n0 and
for all φ ∈ X we have

|DJ(vn)φ| ≤ ε ∥φ∥X . (6.2)

First we prove that {vn}n is bounded in X. Passing to a subsequence of {vn}n (which for
simplicity we denote the same), arguing by contradiction, we assume that ∥vn∥X → +∞.

Using (6.1) and (6.2), we estimate

K +
ε

p+ 1
∥vn∥X ≥ J(vn)−

1

p+ 1
DJ(vn)vn.

Assumption pq > 1 and Remark 2.2 (using, for simplicity, a weaker estimate ∥i∥ ≤ 1) yield

K +
ε

p+ 1
∥vn∥X ≥

( q

q + 1
− 1

p+ 1

)
∥vn∥

q+1
q

X +
( λ

p+ 1
− λ

r + 1

)
∥vn∥r+1

X . (6.3)

Since {vn}n is assumed to be unbounded, for all n sufficiently large ∥vn∥X > 0. Therefore we can

divide both sides of (6.3) by ∥vn∥
q+1
q

X , which gives

K

∥vn∥
q+1
q

X

+
ε

(p+ 1) ∥vn∥1/qX

≥
( q

q + 1
− 1

p+ 1

)
+
( λ

p+ 1
− λ

r + 1

) 1

∥vn∥
1
q−r

X

. (6.4)

From (1.3), we know that all the powers of ∥vn∥X in (6.4) are positive and

q

q + 1
− 1

p+ 1
> 0.

Thus, taking the limit as n → ∞ in (6.4) yields a contradiction and hence {vn}n is bounded.



10 O. AGUDELO, G. HOLUBOVÁ, M. KUDLÁČ EJDE-2025/92

We pass to a subsequence, denoted for simplicity as {vn}n. From Eberlain-Smulyan’s Theorem
(see [19, Th. 2.1.25, p. 67], and from Lemma 2.3, and since X is compactly embedded into C[0, 1],
there exists v ∈ X such that

vn ⇀ v in X, (6.5)

vn → v in L
q+1
q (0, 1). (6.6)

We now prove that {vn}n converges strongly in X. Since (6.5) holds and X is a uniformly
convex space, it follows from [7, Prop. 3.32, p. 78], that it suffices to show that ∥vn∥X → ∥v∥X .
We denote

εn := ∥DJ(vn)∥X∗ := sup
∥φ∥X=1

|DJ(vn)φ|.

Observe that εn ≥ 0 and limn→∞ εn = 0. Choosing φ := vn − v for any n ∈ N, (6.2) reads as
|DJ(vn)φ|

=

∫ 1

0

|v′′n|
1
q−1v′′n(v

′′
n − v′′)dx−

∫ 1

0

(
λ(vn)

r−1
+ vn(vn − v) + (vn)

p−1
+ vn(vn − v)

)
dx

≤ εn∥φ∥X .

Using the triangle inequality, we estimate∣∣ ∫ 1

0

|v′′n|
1
q−1v′′n(v

′′
n − v′′) dx

∣∣− ∣∣ ∫ 1

0

λ(vn)
r−1
+ vn(vn − v) dx

∣∣− ∣∣ ∫ 1

0

(vn)
p−1
+ vn(vn − v) dx

∣∣
≤ εn∥vn − v∥X .

(6.7)

Using that (vn)
r
+ ≤ |vn|r and the Hölder inequality for the conjugate exponents q+1

q and q + 1

yields ∣∣ ∫ 1

0

λ(vn)
r−1
+ vn(vn − v) dx

∣∣ ≤ |λ|∥vn∥rLr(q+1)(0,1)∥vn − v∥L(q+1)/q(0,1).

Since {vn}n is a bounded sequence in X, from Remark 2.2, and (6.6), we obtain

lim
n→∞

∫ 1

0

λ(vn)
r−1
+ vn(vn − v) dx = 0. (6.8)

Arguing in a similar manner yields

lim
n→∞

∫ 1

0

(vn)
p−1
+ vn(vn − v) dx = 0. (6.9)

We recall that {vn}n is bounded in X and εn → 0. Therefore, expressions (6.7), (6.8), and (6.9)
imply

lim
n→∞

∫ 1

0

|v′′n|
1
q−1v′′n(v

′′
n − v′′) dx = 0. (6.10)

In addition, the weak convergence stated in (6.5) yields

lim
n→∞

∫ 1

0

|v′′|
1
q−1v′′(v′′n − v′′) dx = 0 (6.11)

and hence subtracting the terms in (6.10) and (6.11) and using the Hölder inequality yields

0 = lim
n→∞

(
∥vn∥

q+1
q

X −
∫ 1

0

|v′′n|
1
q−1v′′nv

′′ dx−
∫ 1

0

|v′′|
1
q−1v′′v′′n dx+ ∥v∥

q+1
q

X

)
≥ lim

n→∞

(
∥vn∥

q+1
q

X − ∥vn∥
1
q

X∥v∥X − ∥v∥
1
q

X∥vn∥X + ∥v∥
q+1
q

X

)
= lim

n→∞

(
∥vn∥1/qX − ∥v∥1/qX

)(
∥vn∥X − ∥v∥X

)
.

Since the function x 7→ x
1
q is strictly increasing,

0 ≥ lim
n→∞

(
∥vn∥1/qX − ∥v∥1/qX

)
(∥vn∥X − ∥v∥X) ≥ 0;



EJDE-2025/92 ODES WITH CONCAVE-CONVEX NONLINEARITIES 11

thus, necessarily,

∥vn∥X → ∥v∥X . (6.12)

Assertions (6.5) and (6.12) prove the statement. □

7. Proof of Theorem 1.1 and range of λ

We continue the theoretical discussion with the proof of our main result.

Proof of Theorem 1.1. Let λ0 be as in Lemma 4.1 and let us consider λ ∈ (0, λ0). Then it follows
from Propositions 4.3 and 5.2 that there exist two weak nontrivial solutions v1, v2 ∈ X of (3.3)–
(3.4). Moreover, since J(v1) > 0 > J(v2), the weak solutions are necessarily distinct. Using
Proposition 3.1, we verify that v1, v2 are classical solutions of (3.3)–(3.4). Expression (3.2) yields
the corresponding nontrivial smooth functions u1, u2 such that the pairs of functions (u1, v1) and
(u2, v2) solve (3.1). As it was described at the beginning of Section 3, all these functions are
necessarily nonnegative. In conclusion, for λ ∈ (0, λ0), the pairs (u1, v1) and (u2, v2) represent two
distinct nontrivial nonnegative classical solutions to (1.1)-(1.2). This concludes the proof. □

Theorem 1.1 states the existence of at least two distinct solutions of (1.1) for λ ∈ (0, λ0). This
means that in the bifurcation diagram, we find at least two branches of solutions for small positive
values of λ. A natural question concerns with the maximum value of λ0 for which Theorem 1.1 is
still valid.

As we will see in Section 9, when p = 3, q = 1.5 and r = 3−1, the numerical illustration predicts
the existence of the upper and lower branches for λ ∈ (0, λbif), where λbif ≈ 49. For the same
values of the parameters p, q, and r, relation (4.4) provides us with

λ0 ≈ 2.21 ≪ λbif .

This means that the theoretical results herein describe the system only in a narrow interval for λ.
In part, this non-optimality for the range of λ is because the constants in the embeddings stated
in Lemma 2.1 and Remark 2.2 are also not optimal.

To extend the range of values of λ for which Theorem 1.1 is still valid, we track up the energy
estimates in the proof of Lemma 4.1. At this point it is reasonable to use the optimal constant
for the embedding X ↪→ Ls(0, 1) (see inequalities in (2.5)) in Remark 2.2. For instance, similarly
to the well-known result for the Laplace operator, it is readily verified that the infimum

C−1
emb := inf

u∈X, u ̸=0

∥u∥X
∥u∥L(q+1)/q(0,1)

(7.1)

is a positive number and that it is attained. It can also be shown that C−1
emb corresponds to the

principal eigenvalue of the problem

d2

dx2

(
|u′′|γ−2 u′′) = λ|u|γ−2u, x ∈ (0, 1),

u(0) = u(1) = u′′(0) = u′′(1) = 0.

We do not know the exact value of this principal eigenvalue, however, in [5] it is proved that

Cemb ≤ Kemb :=
( 12 )

2

2
min

{( √
πΓ(γ)

Γ(γ + 1
2 )

− 1

γ

)1/γ

,
( √

πΓ(γ′)

Γ(γ′ + 1
2 )

− 1

γ′

)1/γ′}
, (7.2)

where γ = q+1
q , γ′ := γ

γ−1 and Γ(z) :=
∫ +∞
0

tz−1e−tdt.

Let us now improve the estimate on λ0. Since qr < 1, for any v ∈ X

∥v+∥Lr+1(0,1) ≤ ∥v∥Lr+1(0,1) ≤ ∥v∥L(q+1)/q(0,1)

and using (7.1) and (7.2), we obtain that

∥v+∥Lr+1(0,1) ≤ Kemb ∥v∥X . (7.3)



12 O. AGUDELO, G. HOLUBOVÁ, M. KUDLÁČ EJDE-2025/92

From (3.8) we estimate

J(v) ≥ q

q + 1
∥v∥

q+1
q

X −
λKr+1

emb

r + 1
∥v∥r+1

X − 1

2p+1(p+ 1)
∥v∥p+1

X

≥ ∥v∥r+1
X

( q

q + 1
∥v∥

1
q−r

X −
λKr+1

emb

r + 1
− 1

2p+1(p+ 1)
∥v∥p−r

X

)
.

(7.4)

If we replace (4.2) by (7.4) in the proof of Lemma 4.1, we obtain (for the definition of T , see
the proof of Lemma 4.1)

λ0 =
r + 1

Kr+1
emb

( q

q + 1
T

1
q−r − 1

2p+1(p+ 1)
T p−r

)
. (7.5)

When p = 3, q = 1.5, and r = 3−1, we approximate this updated value of λ0 as λ0 ≈ 16.02. Even
though the previous procedure improves the value of λ0, the estimate is still far from optimal.

One of the reasons is that the optimality was used only for the embedding X ↪→ L
q+1
q (0, 1) to

treat the term with Lr+1-norm. The term with Lp+1-norm was treated as before. Another reason
is that solutions predicted in Proposition 4.3 have positive energy, but numerical simulations in
Section 9 indicate there exist solutions of Mountain pass type with nonpositive energy.

8. Symmetry of solutions

Our next result states symmetry properties of solutions of (1.1)-(1.2).

Proposition 8.1. Let λ > 0 and let (u, v) be a classical solution of system (1.1)-(1.2) with u, v > 0
in (0, 1). Then,

(i) u and v are symmetric with respect to the vertical line x = 1/2;
(ii) u

(
1/2

)
= ∥u∥L∞(0,1) and v

(
1/2

)
= ∥v∥L∞(0,1);

(iii) x = 1/2 is the only critical point of u and v.

The proof of Proposition 8.1 follows the (by now) standard method of moving planes (see
[6, 22] and [16, Section 9]). Nonetheless, we present a self-contained proof, adapted to the one
dimensional case treated in this work.

We remark in advance that we only require the monotonicity of the nonlinearities and hence
in this part it is enough to assume 0 < r, p, q < +∞ and λ to be nonnegative. Since we deal with
the one dimensional case, the proof proceeds with an ad-hoc procedure described by a series of
lemmas. The core idea is based upon the classic method of moving planes (see [17, 22]).

Lemma 8.2. Let −∞ ≤ a < b ≤ +∞ and w ∈ C2(a, b) be such that w′′ does not vanish in (a, b),
then w has at most one critical point in (a, b).

The proof of the above lemma follows directly by a direct application of the Mean Value Theorem
to w′.

Lemma 8.3. Let −∞ < a < b < +∞ and w ∈ C2(a, b) ∩ C1[a, b]. Then

(i) w′′ ≤ 0 in (a, b), w(a) = 0 and w(b) ≥ 0 implies that either w ≡ 0 in [a, b] or

w > 0 in (a, b) and w′(a) > 0.

(ii) w′′ ≥ 0 in (a, b), w(a) = 0 and w(b) ≤ 0 implies that either w ≡ 0 in [a, b] or

w < 0 in (a, b) and w′(a) < 0.

Proof. We only prove (i), since (ii) is analogous. First, let x0 ∈ (a, b] be arbitrary, but fixed.

Consider z(x) := w(x) − w(x0)
x0−a (x − a) for x ∈ [a, x0]. Observe that z′′ = w′′ ≤ 0 in (a, x0),

z(a) = w(a) = 0 and z(x0) = 0. The concavity of z implies that z ≥ 0 in (a, x0). Thus,

w(x) ≥ w(x0)
x0−a (x− a) in [a, x0]. In particular, from the definition of derivative, w′(a) ≥ w(x0)

x0−a .

Proceeding similarly one can show that w(x) ≥ w(b) − w(b)−w(x0)
b−x0

(b − x) ≥ 0 in [x0, b] and

w′(b) ≤ w(b)−w(x0)
b−x0

≤ 0. Taking x0 = b, we conclude that w ≥ 0 in [a, b]. Now, let x0 ∈ [a, b]

such that w(x0) = ∥w∥L∞(a,b). If w(x0) = 0, then w ≡ 0 in [a, b]. Otherwise, x0 ∈ (a, b] and



EJDE-2025/92 ODES WITH CONCAVE-CONVEX NONLINEARITIES 13

w(x0) > 0. The above developments imply that w(x) ≥ w(x0)
x0−a (x− a) > 0 in (a, x0), w

′(a) > w(x0)
x0−a

and w(x) ≥ w(b)− w(b)−w(x0)
b−x0

(b− x) > 0 in [x0, b]. This completes the proof. □

Remark 8.4. From the proof of Lemma 8.3, in the case that w′′ ≤ 0 in (a, b), w(a) = w(b) = 0
and ∥w∥L∞(a,b) > 0, we have w′(a) > 0 > w′(b).

Next, consider the space

X0 := C2(0, 1) ∩
{
w ∈ C1[0, 1] : w(0) = w(1) = 0

}
. (8.1)

For w ∈ X0 with w′(0) > 0 and w′(1) < 0, we set

Aw := sup
{
a ∈ (0, 1) : w′ > 0 in (0, a)

}
Bw := inf

{
b ∈ (0, 1) : w′ < 0 in (b, 1)

}
.

Observe that Aw, Bw are well defined and Aw ∈ (0, 1] and Bw ∈ [0, 1). Since w ∈ X0, Aw, Bw ∈
(0, 1).

Lemma 8.5. Let w ∈ X0 be as above. Then Aw and Bw are critical points of w.

Proof. By the approximation property of the supremum and the infimum, we find that w′(Aw) ≥ 0
and w′(Bw) ≤ 0. If w′(Aw) > 0, the continuity of w′ yields that for some δ > 0 small, w′ > 0 in an
interval of the form (0, Aw+δ), thus violating the definition of Aw. We conclude that w′(Aw) = 0.
Similarly, we prove that w′(Bw) = 0. □

Proof of Proposition 8.1. Let u, v ∈ C2(0, 1) ∩ C1[0, 1] with u, v > 0 in (0, 1) and such that the
pair (u, v) is a classical solution of (1.1)-(1.2). Let δ ∈ ( 12 , 1) be arbitrary, but fixed. Observe that
0 < 2δ − 1 < 1. We write

uδ(x) := u(2δ − x) and vδ(x) := v(2δ − x) for x ∈ [2δ − 1, 1].

Notice that uδ and vδ are the reflections of u and v respectively, with respect to the vertical line
x = δ. Also, uδ and vδ solve the system (1.1) in (2δ − 1, 1).

Now we define

wδ(x) := uδ(x)− u(x) and zδ(x) := vδ(x)− v(x) for x ∈ [2δ − 1, 1].

In view of Lemma 8.3 and Remark 8.4, for δ ∈ ( 12 , 1) close enough to 1, u, v are strictly decreasing
in (2δ − 1, 1). Consequently, wδ, zδ > 0 in (δ, 1]. We set

δ∗ := inf
{
δ ∈ (1/2, 1) : zδ > 0 in (δ, 1)

}
.

We claim that

δ∗ = inf
{
δ ∈ (1/2, 1) : wδ > 0 in (δ, 1)

}
.

To prove the claim, let δ ∈ ( 12 , 1) be such that zδ ≥ 0 in [δ, 1]. Using (1.1),

−w′′
δ = λ(vrδ − vr) + vpδ − vp ≥ 0 in (δ, 1).

Since wδ(δ) = 0 and wδ(1) = u(2δ − 1) > 0, Lemma 8.3 implies that wδ > 0 in (δ, 1). This proves
that wδ > 0 in (δ, 1), whenever zδ ≥ 0 in (δ, 1).

Proceeding similarly, zδ > 0 in (δ, 1), whenever wδ ≥ 0 in (δ, 1). Since δ ∈ ( 12 , 1) is arbitrary,
the previous discussion proves the claim.

Now, from the definition of δ∗,

wδ∗(δ∗) = zδ∗(δ∗) = 0.

We prove next that δ∗ = 1
2 . Assume by contradiction that δ∗ > 1

2 and notice that wδ∗ , zδ∗ ≥ 0 in
(δ∗, 1).

Also, since wδ∗(1) = u(2δ∗ − 1) > 0 and zδ∗(1) = v(2δ∗ − 1) > 0, Lemma 8.3 yields that
wδ∗ , zδ∗ > 0 in (δ∗, 1). The continuity of the family of functions {wδ}δ and {zδ}δ in C1[0, 1] with

respect to the parameter δ, allows us to find δ̂ ∈ (1/2, δ∗) such that wδ̂, zδ̂ > 0 in (δ̂, 1). This
contradicts the definition of δ∗ and proves that δ∗ = 1/2.



14 O. AGUDELO, G. HOLUBOVÁ, M. KUDLÁČ EJDE-2025/92

Since w 1
2
(x) = u(1 − x) − u(x) ≥ 0 and z 1

2
(x) = v(1 − x) − v(x) ≥ 0 for every x ∈ [ 12 , 1], we

find that u(1 − x) ≥ u(x), v(1 − x) ≥ v(x) for any x ∈ [ 12 , 1]. A similar argument shows that

u(1− x) ≥ u(x), v(1− x) ≥ v(x) for any x ∈ [0, 1
2 ] and consequently for x ∈ [0, 1].

Now notice that u(1 − x), v(1 − x) also solve (1.1)-(1.2). Since (u 1
2
) 1

2
= u and (v 1

2
) 1

2
= v, we

may argue as above to find that u(1− x) ≤ u(x), v(1− x) ≤ v(x) for any x ∈ [0, 1]. Thus, u, v are
symmetric with respect to x = 1

2 .
Now, let xu, xv ∈ [0, 1] be such that

u(xu) = ∥u∥L∞(0,1) and v(xv) = ∥v∥L∞(0,1).

Since u(xu) = max{u(x) : x ∈ [0, 1]} and u′(0) > 0 > u′(1) ̸= 0, xu ∈ (0, 1) and consequently
u′(xu) = 0. Similarly, xv ∈ (0, 1) and v′(xv) = 0.

From Lemmas 8.2 and 8.5 and the symmetry of u and v, Au = Bu = xu = 1
2 and Av = Bv =

xv = 1
2 and completes the proof. □

9. Numerical illustration

To support the theoretical results in this paper and to obtain wider intuition about behavior
of our system, this section discusses the numerical strategy for obtaining visual description of our
system.

Let us begin with the strategy for the implementation. First, we transform the system (1.1)
into the associated first order system (of four equations)

u′(x) = w(x),

v′(x) = z(x),

w′(x) = −λ(v+(x))
r − (v+(x))

p,

z′(x) = −|u(x)|q−1 u(x),

(9.1)

and consider this system with the initial conditions

u(0) = 0, v(0) = 0, w(0) = du0, z(0) = dv0, (9.2)

where du0 and dv0 are considered as free parameters. Let us assume that the initial boundary value
problem (9.1)-(9.2) is such that for (du0, dv0) in a rectangle R ⊂ (0,+∞) × (0,+∞), existence,
uniqueness and continuity of solutions with respect to the parameters (du0, dv0) hold in [0, 1].
Then given (du0, dv0) ∈ R, the functions u, v, w and z are continuous in [0, 1] and we write

u(x) = u(x, du0, dv0) and v(x) = v(x, du0, dv0).

We are therefore interested in the zeroes of the continuous mapping Φ : R → R2 defined by

Φ(du0, dv0) = (u(1), v(1)). (9.3)

Analyzing the sign changes of the components of Φ, with respect to the values of du0 and dv0, we
locate numerically two solutions for a wide range of values λ. As mentioned in the Introduction
such implementation is motivated by a combination of the standard Shooting method and the
heuristics of the Poincaré-Miranda Theorem.

The strategy is as follows.

• Choose du0 and dv0 appropriately, i.e., so that existence and uniqueness in [0, 1] and
continuity with respect to initial conditions holds.

• Compute the corresponding solutions (u, v, w, z) of (9.1)-(9.2) in the interval [0, 1].
• Extract the values u(1) and v(1). Observe that for certain values of du0 and dv0 the
corresponding pair of functions (u, v) is a solution to (3.1) provided u(1) = v(1) = 0.

• Given a tolerance ϵ > 0, the numerical calculations of u and v yield an admissible numerical
solution to (3.1) provided |u(1)|, |v(1)| < ϵ.

• It is more efficient to track simultaneously the different regions where either |u(1)| ≥ ϵ or
|v(1)| ≥ ϵ. This can be interpreted as tracking the sign-changes of the values u(1), v(1) as
the parameters du0 and dv0 vary.



EJDE-2025/92 ODES WITH CONCAVE-CONVEX NONLINEARITIES 15

A reader can imagine the process as if two football players kick simultaneously two balls from
the ground (zero height) on the left border of the field, each of the kicks performed with a given
initial slope. The trajectory of the kicked balls are ruled by the equations in (3.1). The football
players aim at a bin located on the right border of the field. To hit the bin with a ball it is
necessary to choose the initial slope of the kick so that the ball descends to a zero height exactly
on the right border of the field. Since each of the football players may kick with different strength,
the required slopes for the football players need not be the same.

The implementation initially uses a “coarse” grid of different parameters du0, dv0 to roughly
locate the sign-changes. Then, an adaptive strategy is employed in order to improve the grid
density. This strategy iterates the values of the parameters du0 and dv0, computes again the
corresponding solutions of (9.1) and tracks the regions of the corresponding changes of sign.

We implemented the idea of the experiments as a set of scripts in Matlab, see [26] for the
implementation. The scripts are optimized to provide the results in reasonable time and accuracy.

<numerical.m>:
Define algorithm settings and values of the parameters p, q, r.
Define the initial range for du0 and dv0.
Define the range for the parameter λ.
Iterate through the range for λ and do the following:
<for cycle>:

Run shooting.m with desired settings.
<shooting.m>:

Represent du0-dv0 plane by coarse (∆ = 0.1) and dense (∆ = 0.005) grids.
Iterate vertices of the coarse grid and run shootandsolve.m.
<shootandsolve.m>:

Run ode45.m (Runge-Kutta method for ODEs from Matlab library) for
(9.1)-(9.2) and any given initial condition
(u(0), w(0), v(0), z(0)) = (0, du0, 0, dv0).
Compute the solution and calculate the residues u(1) and v(1).
Return the residues as the return value.

</shootandsolve.m>:
Based on the residues for a pair (du0, dv0), assign the pair a color using the
following scheme:

u(1) > 0, v(1) > 0 − green;

u(1) > 0, v(1) < 0 − yellow;

u(1) < 0, v(1) > 0 − blue;

u(1) < 0, v(1) < 0 − red.

(9.4)

Plot the colors in a du0-dv0 diagram and save them into a variable for the coarse
grid.
The pair (du0, dv0) such that (u, v) is also a solution of (3.1) can be located
exactly at the point where all colors meet, in other words, where both residues
are zero.
Iterate through the vertices of the coarse grid and choose only the points which
have a neighbor of a different color (we target edges between two colors).
In the dense grid, proceed only with the vertices corresponding to the chosen
points in the coarse grid.
Run /shootandsolve.m for the vertices in the dense grid.
Plot the colors in the du0-dv0 diagram and save them into a variable for the
dense grid.
Explore neighbourhood of the points in the dense grid and check whether all
colors are present in the neighbourhood (if so, assume there is a solution in the
neighbourhood).



16 O. AGUDELO, G. HOLUBOVÁ, M. KUDLÁČ EJDE-2025/92

Approximate the values of (du0, dv0) corresponding to the solution – denote
them by (solU, solV ) – and mark them in the graph with a black circle.
Return (solU, solV ) – or (Inf, Inf) if no solution was found.

</shooting.m>:
If the return value is (Inf, Inf) (solution not found), exit the “for” cycle.
Run showsolution.m, pass (solU, solV ) as a parameter.

<showsolution.m>:
Compute the solution using (solU, solV ) and ode45.m.
Plot ∥v∥L∞(0,1) vs λ, where v is the component v found above, i.e., the numerical
solution of (3.3).
Return L∞-norm of the plotted function.

</showsolution.m>:
Save the value of λ and the corresponding ∥v∥L∞(0,1) of the solution v into a text file.
Considering the development of the results for the previous values of λ, automatically
adapt the ranges of du0 and dv0 for the next iteration (for the optimization, it is
necessary to use the narrowest range of the parameters as possible).

</for cycle>:
Load results for the whole range of λ from the text file.
Vizualize dependancy of the L∞-norm of the solution on the parameter λ – plot the

bifurcation diagram.
</numerical.m>:

The grid is rectangled and uniform, and the distance (in either direction du0 and dv0) between
two neighboring vertices is ∆. The results are computed only in the vertices of the grid. For
instance, consider the rectangle [a, b] × [c, d] ⊂ (0,∞) × (0,∞) for selecting the pairs (du0, dv0).
The corresponding grid of vertices reads as{

(a+ i∆, b+ j∆) ∈ [a, b]× [c, d] : i = 0, 1, . . . ,
b− a

∆
, j = 0, 1, . . . ,

d− c

∆

}
.

As an example, if b− a, d− c,∆−1 ∈ N, then for any unit square in du0-dv0 plane, the coarse grid
contains 102 vertices and the dense grid contains 2002 vertices.

The use of a coarse grid first and then a denser one in the implementation saves significant
amount of time and memory. In the graph, the optimized script skips parts of the domain where
the results cannot be located (from the coarse grid’s point of view), thus the script leaves blank
rectangles in the graph.

Next, we discuss our numerical findings in more detail. Recall that we have set p = 3, q = 1.5,
and r = 1/3. Nevertheless, for small positive values of λ, the numerical experiments anticipate
existence of a nontrivial solution as it can be seen in Figure 2.

Figure 2b shows a narrow red protrusion coming from the red area at the bottom-left corner of
the diagram. At the point where the protrusion touches the green area, all four colors connect at
the pair (du0, dv0) corresponding to a nontrivial numerical solution.

As the value of λ increases, the red protrusion is visible for higher and higher values of du0 and
dv0, as shown in Figure 3a.

If we fix λ = 10, besides the solution illustrated in Figure 3a with du0 ≈ 1 and dv0 ≈ 0.03,
the experiments found another solution for du0 ≈ 44, dv0 ≈ 16.5, as illustrated in Figure 3b.
Also, Figures 3c, 3d show that the L∞-norm of the solution for lower initial slopes du0, dv0
(corresponding to the diagram in Figure 3a) is lower than the L∞-norm of the other solution. The
numerical experiments yield a similar scenario when λ = 1.

The numerical results discussed above strongly suggest the existence of two branches in the
bifurcation diagram; the lower branch (closer to the trivial solution – in the view of the L∞-
norm) and the upper branch (farther from the trivial solution). Based on the presented results,
the branches are also getting closer to each other with increasing λ, presumably colliding when
λ = λbif .

With this assumption, we let the script explore the range 1 ≤ λ ≤ 50. The results confirmed
the assumptions that the branches meet at some point λbif ≈ 49. The results also confirmed



EJDE-2025/92 ODES WITH CONCAVE-CONVEX NONLINEARITIES 17

(a) Diagram for λ = 1 (coarse grid) (b) Diagram for λ = 1 (dense grid). The
pair of parameters corresponding to a so-
lution is marked by a black dot

Figure 2. du0-dv0 diagram for λ = 1 and both coarse and dense grid

(a) du0-dv0 diagram for
λ = 10, du0 ≈ 1, dv0 ≈ 0.03

(b) du0-dv0 diagram for
λ = 10, du0 ≈ 44, dv0 ≈ 16.5

(c) Lower branch solution for
λ = 10 corresponding to Figure 3a

(d) Upper branch solution for
λ = 10 corresponding to Figure 3b

Figure 3. Numerical experiments for λ = 10

(numerically) that beyond λbif , no solutions could be found. This behavior is summarized in the
bifurcation diagram in Figure 1. Further illustrations can be found in Figures 4 and 5.

Acknowledgements. The authors were supported by the Grant CR 22-18261S from the Grant
Agency of the Czech Republic.



18 O. AGUDELO, G. HOLUBOVÁ, M. KUDLÁČ EJDE-2025/92

(a) Solution v5 in du0-dv0 dia-
gram

(b) Solution v20 in du0-dv0 dia-
gram

(c) Solution v40 in du0-dv0 dia-
gram

Figure 4. Solutions from upper branch of the bifurcation diagram shown in du0-
dv0 diagram

(a) Solution v5 in du0-dv0 dia-
gram

(b) Solution v20 in du0-dv0 dia-
gram

(c) Solution v40 in du0-dv0 dia-
gram

Figure 5. Solutions from lower branch of the bifurcation diagram shown in du0-
dv0 diagram

References

[1] Robert A. Adams, John J. F. Fournier; Sobolev spaces, second ed., Elsevier, 2003.

[2] Oscar Agudelo, Bernhard Ruf, Carlos Veléz; On a hamiltonian elliptic system with concave and convex non-
linearities, Discrete and Continuous Dynamical Systems-S 16 (2023), no. 11, 2902–2918.

[3] Oscar Agudelo, Bernhard Ruf, Carlos Velez; Multiplicity results for a subcritical hamiltonian system with

concave-convex nonlinearities, arXiv preprint arXiv:2412.10812 (2024), Accepted in Calc. Var.
[4] A. Ambrosetti, H. Brezis, G. Cerami; Combined effects of concave and convex nonlinearities in some elliptic

problems, J. Funct. Anal. 122 (1994), no. 2, 519–543.
[5] J. Benedikt; Estimates of the principal eigenvalue of the p-laplacian and the p-biharmonic operator, Mathe-

matica Bohemica 140 (2015), no. 2, 215–222.
[6] Henri Berestycki, Louis Nirenberg; Monotonicity, symmetry and antisymmetry of solutions of semilinear

elliptic equations, Journal of Geometry and Physics 5 (1988), no. 2, 237–275.

[7] Haim Brezis; Functional analysis, sobolev spaces and partial differential equations, Springer, 2011.

[8] A. Castro, A. C. Lazer; On periodic solutions of weakly coupled systems of differential equations, Bollettino
dell’Unione Matematica Italiana, Serie B 18 (1981), no. 3, 733–742 (Italian).

[9] Alfonso Castro, Jorge Cossio, Sigifredo Herrón, Carlos Vélez; Shooting from singularity to singularity and a
quasilinear p-laplace-beltrami equation with indefinite weight, Discrete and Continuous Dynamical Systems-S
17 (2024), no. 5&6, 2186–2207.

[10] Alfonso Castro, Alexandra Kurepa; Radial solutions to a dirichlet problem involving critical exponents when

n = 6, Transactions of the American Mathematical Society 348 (1996), no. 2, 781–798.
[11] Alfonso Castro, Vı́ctor Padrón; Classification of radial solutions arising in the study of thermal structures

with thermal equilibrium or no flux at the boundary, vol. 208, American Mathematical Society, 2010.
[12] Alfonso Castro, R Shivaji; Multiple solutions for a dirichlet problem with jumping nonlinearities. ii, Journal

of mathematical analysis and applications 133 (1988), no. 2, 509–528.

[13] S Chandrasekhar; An introduction to the study of stellar structure, The University of Chicago Press, 1939.
[14] Philippe Clément, Patricio Felmer, Enzo Mitidieri; Homoclinic orbits for a class of infinite dimensional hamil-

tonian systems, Annali della Scuola Normale Superiore di Pisa-Classe di Scienze 24 (1997), no. 2, 367–393.
[15] Francesca Dalbono, PJ McKenna; Multiplicity results for a class of asymmetric weakly coupled systems of

second-order ordinary differential equations, Boundary Value Problems 2005 (2005), 1–23.



EJDE-2025/92 ODES WITH CONCAVE-CONVEX NONLINEARITIES 19

[16] Djairo G. De Figueiredo; Semilinear elliptic systems: existence, multiplicity, symmetry of solutions, Handbook

of differential equations: stationary partial differential equations 5 (2008), 1–48.

[17] Djairo G. De Figueredo; Monotonicity and symmetry of solutions of elliptic systems in general domains,
Nonlinear Differential Equations and Applications NoDEA 1 (1994), no. 2, 119–123.

[18] E. M. dos Santos; On a fourth-order quasilinear elliptic equation of concave-convex type, NoDEA Nonlinear

Differential Equations Appl. (2009).
[19] P. Drábek, J. Milota; Methods of nonlinear analysis: applications to differential equations, Springer Science

& Business Media, 2007.

[20] Lawrence C. Evans; Partial differential equations, vol. 19, American Mathematical Soc., 2010.
[21] Svatopluk Fućık, Alois Kufner; Nonlinear differential equations, Elsevier, 1980.

[22] Basilis Gidas, Wei-Ming Ni, Louis Nirenberg; Symmetry and related properties via the maximum principle,

Communications in mathematical physics, 68 (1979), no. 3, 209–243.
[23] Josephus Hulshof, Robertus van der Vorst; Asymptotic behaviour of ground states, Proceedings of the American

Mathematical Society 124 (1996), no. 8, 2423–2431.
[24] R. Kippenhaln, A. Weigert; Stellar structure and evolution, Springer-Verlag, 1990.

[25] Systems with concave-convex nonlinearities: existence and multiplicity of solutions, Diploma thesis,

Západoeská univerzita Plzni, 2023.
[26] M. Kudláč; Poincaremirandashooting: Scripts for finding solutions of a system of two odes with boundary

conditions, https://github.com/Arkandrus/PoincareMirandaShooting, 2025.

[27] Mike R. Osborne; On shooting methods for boundary value problems, Journal of mathematical analysis and
applications 27 (1969), no. 2, 417–433.

[28] Bernhard Ruf; Superlinear elliptic equations and systems, Handbook of differential equations: stationary

partial differential equations 5 (2008), 211–276.
[29] Michael Struwe; Variational methods, fourth ed., vol. 34, Springer, 2008.

[30] Vladimir A. Zorich; Mathematical Analysis 1, Springer, 2004.

Oscar Agudelo
Department of Mathematics and NTIS, Faculty of Applied Sciences, University of West Bohemia in

Pilsen, Univerzitńı 8, 301 00 Plzeň, Czech Republic

Email address: oiagudel@kma.zcu.cz, ORCID: 0000-0002-2588-9999

Gabriela Holubová

Department of Mathematics and NTIS, Faculty of Applied Sciences, University of West Bohemia in
Pilsen, Univerzitńı 8, 301 00 Plzeň, Czech Republic

Email address: gabriela@kma.zcu.cz, ORCID: 0000-0003-1127-3381

Martin Kudláč
Department of Mathematics and NTIS, Faculty of Applied Sciences, University of West Bohemia in

Pilsen, Univerzitńı 8, 301 00 Plzeň, Czech Republic
Email address: kudlacm@kma.zcu.cz, ORCID: 009-0007-1749-3314


	1. Introduction
	2. Preliminaries
	3. Functional analytic setting
	4. Mountain pass solution
	5. Local minimum solution
	6. Palais-Smale condition
	7. Proof of Theorem 1.1 and range of 
	8. Symmetry of solutions
	9. Numerical illustration
	Acknowledgements

	References