Electronic Journal of Differential Equations, Vol. 2025 (2025), No. 43, pp. 1–21. ISSN: 1072-6691. URL: https://ejde.math.txstate.edu, https://ejde.math.unt.edu DOI: 10.58997/ejde.2025.43 EXISTENCE OF MAXIMAL AND MINIMAL WEAK SOLUTIONS AND FINITE DIFFERENCE APPROXIMATIONS FOR ELLIPTIC SYSTEMS WITH NONLINEAR BOUNDARY CONDITIONS SHALMALI BANDYOPADHYAY, THOMAS LEWIS, NSOKI MAVINGA Abstract. We establish the existence of maximal and minimal weak solutions between or- dered pairs of weak sub- and super-solutions for a coupled system of elliptic equations with quasimonotone nonlinearities on the boundary. We also formulate a finite difference method to approximate the solutions and establish the existence of maximal and minimal approximations between ordered pairs of discrete sub- and super-solutions. Monotone iterations are formulated for constructing the maximal and minimal solutions when the nonlinearity is monotone. Nu- merical simulations are used to explore existence, nonexistence, uniqueness and non-uniqueness properties of positive solutions. When the nonlinearities do not satisfy the monotonicity condi- tion, we prove the existence of weak maximal and minimal solutions using Zorn’s lemma and a version of Kato’s inequality up to the boundary. 1. Introduction We consider the coupled system of elliptic equations with nonlinear boundary conditions −∆ui + ui = 0 in Ω; ∂ui ∂η = fi(x, u1, u2) on ∂Ω, i = 1, 2, (1.1) where Ω ⊂ RN is a bounded domain with Lipschitz (C0,1) boundary ∂Ω, N ≥ 2, and ∂/∂η := η(x)· ∇ denotes the outer normal derivative on the boundary ∂Ω. For each i = 1, 2, fi : ∂Ω× R2 → R is a Carathéodory function; that is, fi(·, u1, u2) is measurable for all (u1, u2) ∈ R2 and fi(x, ·, ·) is continuous for a.e. x ∈ ∂Ω. Throughout this article we assume that each fi satisfies the quasimonotonicity condition (A1) the functions fi are quasimonotone nondecreasing in the sense that f1(x, u1, u2) is nonde- creasing in u2 for all fixed x ∈ ∂Ω, u1 ∈ R, and f2(x, u1, u2) is nondecreasing in u1 for all fixed x ∈ ∂Ω, u2 ∈ R. In this article, we establish the existence of maximal and minimal weak solutions for (1.1). In particular, we use a monotone iteration method when the fi’s are monotone nondecreasing in both variables u1 and u2, and for the fi’s that are nonmonotone in one of the variables, we utilize the surjectivity of a bounded, pseudomonotone and coercive operator, Zorn’s lemma and a version of Kato’s inequality up to the boundary to obtain the existence of maximal and minimal solutions. To visualize solutions, we utilize numerical methods. Some common numerical approximation techniques for reaction diffusion equations can be found in [19, 20, 23, 24, 25]. We use the finite difference method to approximate solutions inspired by the results for approximating semilinear elliptic problems with Dirichlet boundary conditions in [15]. We formulate a finite difference method for (1.1) and prove admissibility and stability results. We also formulate a methodology 2020 Mathematics Subject Classification. 35J60, 35J67, 65N06, 65N22. Key words and phrases. Weak solutions; quasimonotone; subsolution; supersolution; Zorn’s lemma; finite difference method; Kato’s inequality. ©2025. This work is licensed under a CC BY 4.0 license. Submitted April 17, 2024. Published Aril 23, 2025. 1 2 S. BANDYOPADHYAY, T. LEWIS, N. MAVINGA EJDE-2025/43 for finding maximal and minimal solutions in the monotone case. This methodology is used to generate bifurcation diagrams for several one-dimensional examples. The numerical study complements the analytical results. Elliptic equations are important for studying mathematical models in problems such as chemical reactions, ecology, population dynamics and combustion theory. Extensive studies have been done when the elliptic equation has linear boundary conditions which includes Dirichlet, Neumann, and Robin boundary conditions. However, there are certain scenarios where chemical reactions, the biological bonding, or species interactions may occur in a narrow layer or region near the boundary. In such cases, linear boundary conditions are deficient to describe the mathematical model (see e.g. [5, 10, 13, 16, 22, 21] and the references therein). Therefore, studying elliptic problems with nonlinear boundary conditions has been of much interest over the past few decades (see e.g. [2, 3, 4, 5, 9, 11, 16, 17, 26] and the references therein). In [4] the authors studied the existence of a maximal and a minimal weak solution between ordered pairs of sub and supersolutions for elliptic (scalar) equations with nonlinear boundary conditions for both monotone and nonmonotone nonlinearities. The aim of this paper is to extend the results to the case of elliptic systems with quasimonotone nonlinearities on the boundary, and approximate these solutions using the finite difference method. We shall point out that quasi- monotone elliptic systems have been studied by several authors, we refer to [7, 18, 22]. In [18], the authors proved the existence of a maximal weak solution between ordered pairs of sub and super- solutions for quasimonotone elliptic systems with linear boundary conditions, namely, Dirichlet boundary conditions. Furthermore, they assumed that the nonlinearities are not necessarily dif- ferentiable or even continuous. In such a case the monotone iteration procedure is not applicable, and the main ingredient is the use of pseudomonotone operators theory. As for [7] the author extended the results in [18] to the p-Laplacian case but still with Dirichlet boundary conditions. The author in [22] considered quasimonotone elliptic systems with smooth nonlinear boundary conditions, and they proved the existence of maximal and minimal classical solutions using mono- tone iteration methods. Here, we address the existence of maximal and minimal weak solutions for (1.1) assuming that the nonlinearities are Carathéodory. Throughout this paper H1(Ω) denotes the usual real Sobolev space of functions on Ω; the product Sobolev space H1(Ω) × H1(Ω) will be denoted by (H1(Ω))2 and is endowed with the norm ∥(u1, u2)∥(H1(Ω))2 := ∥u1∥H1(Ω) + ∥u2∥H1(Ω). Moreover, the product space (H1(Ω))2 is reflexive as H1(Ω) is reflexive (see e.g. [12, p.15]). Besides the Sobolev spaces, we use the real Lebesgue space Lq(∂Ω), and the compactness of the trace operator Γ : H1(Ω) → Lq(∂Ω) with Γu = u|∂Ω (see e.g. [1, 8], [6, Thm 2.79], and [12, Chapter 6]), that is, Γ is continuous (compact) if 1 ≤ q ≤ 2(N − 1) N − 2 (1 ≤ q < 2(N − 1) N − 2 ), if N > 2, 1 ≤ q (1 ≤ q) if N = 2. (1.2) To keep the notation simple, we will use the following: U := (u1, u2), U := (u1, u2) and U := (u1, u2). The inequality U ≤ V means ui(x) ≤ vi(x) a.e. x ∈ Ω and a.e. x ∈ ∂Ω for each i = 1, 2. Definition 1.1. We say that a function (u1, u2) ∈ (H1(Ω))2 is a weak solution to (1.1) if (i) For every i = 1, 2, fi(·, u1(·), u2(·)) ∈ Lr(∂Ω) for some r > 1 ifN = 2 and fi(·, u1(·), u2(·)) ∈ L 2(N−1) N (∂Ω) if N > 2 and (ii) ∫ Ω (∇ui∇ψ + uiψ) = ∫ ∂Ω fi(x, u1, u2)ψ for all ψ ∈ H1(Ω). Definition 1.2. We say that a function (u1, u2) ∈ (H1(Ω))2 is a weak supersolution to (1.1) if (i) For every i = 1, 2, fi(·, u1(·), u2(·)) ∈ Lr(∂Ω) for some r > 1 ifN = 2 and fi(·, u1(·), u2(·)) ∈ L 2(N−1) N (∂Ω) if N > 2 and (ii) ∫ Ω (∇ui∇ψ + uiψ) ≥ ∫ ∂Ω fi(x, u1, u2)ψ for all 0 ≤ ψ ∈ H1(Ω). A weak subsolution (u1, u2) is defined by reversing the inequality in (ii) above. Observe that the integrals on the right hand side of (ii) of Definition 1.1 and Definition 1.2 make sense since (i) holds and r = 2(N−1) N is the conjugate of 2(N−1) N−2 when N > 2. EJDE-2025/43 EXISTENCE OF MAXIMAL AND MINIMAL WEAK SOLUTIONS 3 In what follows, we state our main results which include the case where the fi’s are monotone nondecreasing in both variables u1 and u2 (see Theorem 1.3) and the case where the fi’s are nonmonotone in one of the variables (see Theorem 1.4). Theorem 1.3. Assume that there exist a weak subsolution (u1, u2) and a weak supersolution (u1, u2) of (1.1) such that (u1, u2) ≤ (u1, u2) a.e. on Ω, and (A2) there exists k1 > 0 such that for all x ∈ ∂Ω and (s1, s2) ∈ R2 with (u1(x), u2(x)) ≤ (s1, s2) ≤ (u1(x), u2(x)), f1(x, s1, s2) + k1s1 is nondecreasing in s1. (A3) there exists k2 > 0 such that for all x ∈ ∂Ω and (s1, s2) ∈ R2 with (u1(x), u2(x)) ≤ (s1, s2) ≤ (u1(x), u2(x)), f2(x, s1, s2) + k2s2 is nondecreasing in s2. Then, there exist a minimal weak solution (u1,∗, u2,∗) and a maximal weak solution (u∗1, u ∗ 2) to (1.1), that is, if (u1, u2) is any weak solution to (1.1) such that (u1, u2) ≤ (u1, u2) ≤ (u1, u2), then (u1,∗, u2,∗) ≤ (u1, u2) ≤ (u∗1, u ∗ 2). In the next theorem, we remove the monotonicity conditions (A2) and (A3) and obtain the following result. Theorem 1.4. Assume that there exists a pair of ordered weak subsolution (u1, u2) and superso- lution (u1, u2) of (1.1), and that the following condition holds: (A4) there exist K1,K2 ∈ Lr(∂Ω), r > 2(N−1) N , such that |fi(x, s1, s2)| ≤ Ki(x) a.e. x ∈ ∂Ω, whenever ui(x) ≤ si ≤ ui(x), i = 1, 2. Then there exists a weak solution (u1, u2) of (1.1) such that (u1, u2) ≤ (u1, u2) ≤ (u1, u2). More- over, there exist a minimal weak solution (u1,∗, u2,∗) and a maximal weak solution (u∗1, u ∗ 2) to (1.1); that is, for any weak solution (u1, u2) to (1.1) with (u1, u2) ≤ (u1, u2) ≤ (u1, u2), we have (u1,∗, u2,∗) ≤ (u1, u2) ≤ (u∗1, u ∗ 2). In the course of the proofs, we will need the following result on the existence of weak solutions for the single equations case. Proposition 1.5. Consider the nonlinear problem −∆u+ u = 0 in Ω ; ∂u ∂η = f(x, u) on ∂Ω, (1.3) where f : ∂Ω × R → R is a Carathéodory function. Suppose that there exists a pair of a weak subsolution u and a weak supersolution u such that u ≤ u in Ω, and that there exists K ∈ Lr(∂Ω), r > 2(N−1) N , such that |f(x, s)| ≤ K(x) a.e. x ∈ ∂Ω and for all s ∈ R satisfying u(x) ≤ s ≤ u(x). Then (1.3) has at least one weak solution u such that u ≤ u ≤ u. The proof of Proposition 1.5 can be found in [4], which relies on the surjectivity of bounded, coercive pseudomonotone operators, Zorn’s lemma and a version of Kato’s inequality up to the boundary for single equations. This article is organized as follows. In section 2, we prove Theorem 1.3. Section 3 is devoted to the proof of Theorem 1.4. In Section 4, we formulate and analyze a finite difference method to approximate solutions and generate bifurcation diagrams for several one-dimensional cases. Finally, in Section 5, the appendix, we state a version of Kato’s inequality up to the boundary for single equations and employ it to prove that the componentwise maximum of two solutions of (1.1) is a subsolution of (1.1) (see Proposition 5.4) and componentwise minimum of two solutions of (1.1) is a supersolution of (1.1) (see Proposition 5.6), results necessary to prove Theorem 1.4. 2. Proof of Theorem 1.3 We will first construct a monotone operator and then use a corresponding iterative scheme to show the existence of a minimal (maximal) solution using the convergence of a sequence of weak subsolutions (supersolutions). 4 S. BANDYOPADHYAY, T. LEWIS, N. MAVINGA EJDE-2025/43 (I) Construction of the monotone operator. We define the map T : J → (H1(Ω))2 by T (U) =W , where J := {U = (u1, u2) ∈ (H1(Ω))2 : U ≤ U ≤ U} and W = (w1, w2) is the unique weak solution of the decoupled system −∆wi + wi = 0 in Ω; ∂wi ∂η + kwi = fi(x, u1, u2) + kui on ∂Ω, i = 1, 2, (2.1) where k = k1 + k2 ≥ 0. Observe that, for each i, (2.1) is a linear equation for wi. Also notice that fi(x, u1, u2) + kui ∈ L 2(N−1) N (∂Ω). Indeed, fi(,̇u1(·), u2(·)) ∈ L 2(N−1) N (∂Ω) by Definition 1.1 and (1.2) implies ui ∈ L 2(N−1) N (∂Ω) since 2(N−1) N−2 > 2(N−1) N . Hence, by the existence results for linear elliptic equations (see e.g. [14, P.160-162], there exists a unique solution of (2.1). Thus, T is well-defined. Now, we prove that T is monotonically nondecreasing and maps J into itself. Indeed, let U, V ∈ J with U ≤ V . We have that T (U) = W = (w1, w2) and T (V ) = Z = (z1, z2) satisfy the following: −∆wi + wi = 0 in Ω; ∂wi ∂η + kwi = fi(x, u1, u2) + kui on ∂Ω, i = 1, 2, (2.2) and −∆zi + zi = 0 in Ω, ∂zi ∂η + kzi = fi(x, v1, v2) + kvi on ∂Ω, i = 1, 2. (2.3) By (A2) and the fact that (u1, u2) ≤ (v1, v2) and k > 0, we have f1(x, u1, u2) + ku1 ≤ f1(x, v1, u2) + kv1. Applying the quasimonotonicity condition (A1), we conclude f1(x, u1, u2) + ku1 ≤ f1(x, v1, v2)+kv1. Similarly, by (A1) and (A3), we have f2(x, u1, u2)+ku2 ≤ f2(x, v1, v2)+ kv2. Substracting (2.2) and (2.3), we obtain that wi − zi satisfies −∆(wi − zi) + (wi − zi) = 0 in Ω; ∂(wi − zi) ∂η + k(wi − zi) ≤ 0 on ∂Ω. It follows from the comparison principle (see e.g. [8]) that wi ≤ zi. So, T (U) ≤ T (V ). Now, we show that T (J) ⊆ J . With the monotonicity of T , it is sufficient to show that U ≤ T (U) and T (U) ≤ U . Take the subsolution U = (u1, u2). Then T (U) = (w1, w2) satisfies the system −∆wi + wi = 0 in Ω; ∂wi ∂η + kwi = fi(x, u1, u2) + kui on ∂Ω, i = 1, 2. Using the fact that U is a subsolution to (1.1), we obtain that −∆(ui − wi) + (ui − wi) = 0 in Ω; ∂(ui − wi) ∂η + k(ui − wi) ≤ 0 on ∂Ω, i = 1, 2. By the comparison principle, we have ui−wi ≤ 0 in Ω. Hence, U ≤ T (U). Similarly, we can show T (U) ≤ U . Hence, U ≤ T (U) ≤ T (U) ≤ U. (2.4) Thus, T maps J into itself. (II) Construction of minimal and maximal weak solutions. We construct a monotone sequence {Un} = {(u1,n, u2,n)} and a monotone sequence {Wn} = {(w1,n, w2,n)} using the linear iteration process as follows: Un = T (Un−1) with U0 = U, and Wn = T (Wn−1) with W0 = U, EJDE-2025/43 EXISTENCE OF MAXIMAL AND MINIMAL WEAK SOLUTIONS 5 where Un and Wn are weak solutions of −∆ui,n + ui,n = 0 in Ω; ∂ui,n ∂η + kui,n = fi(x, u1,n−1, u2,n−1) + kui,n−1 on ∂Ω, (2.5) and −∆wi,n + wi,n = 0 in Ω; ∂wi,n ∂η + kwi,n = fi(x,w1,n−1, w2,n−1) + kwi,n−1 on ∂Ω for each i = 1, 2. The monotonicity of T and (2.4) imply that U = U0 ≤ U1 ≤ U2 ≤ · · · ≤Wn ≤Wn−1 ≤ · · · ≤W0 = U. We now proceed to show that the sequences Un and Wn are weakly convergent. Since Un = T (Un−1) is a weak solution of (2.5), we have that∫ Ω (∇ui,n∇ψ + ui,nψ) + k ∫ ∂Ω ui,nψ = ∫ ∂Ω ( fi(x, u1,n−1, u1,n−1) + kui,n−1 ) ψ (2.6) for all ψ ∈ H1(Ω). Taking ψ = ui,n for each i = 1, 2, we obtain∫ Ω ( |∇ui,n|2 + u2i,n ) + k ∫ ∂Ω u2i,n = ∫ ∂Ω ( fi(x, u1,n−1, u1,n−1) + kui,n−1 ) ui,n. (2.7) Observe that ∥fi(x, u1,n−1, u2,n−1) + kui,n−1∥ L 2(N−1) N (∂Ω) ≤ ∥fi(x, u1, u2) + kui∥ L 2(N−1) N (∂Ω) + ∥fi(x, u1, u2) + kui∥ L 2(N−1) N (∂Ω) ≤ C̃ , (2.8) where C̃ is a constant independent of n. Hence, for each i = 1, 2, fi(·, u1,n−1(·), u2,n−1(·)) + kui,n−1(·) ∈ L 2(N−1) N (∂Ω). Then, the right- hand side of (2.7) can be estimated using Hölder’s inequality and the bound (2.8), and so (2.7) can be estimated as follows: ∥ui,n∥2H1(Ω) ≤ ∥ui,n∥2H1(Ω) + k∥ui,n∥2L2(∂Ω) ≤ ∥fi(x, u1,n−1, u2,n−1) + kui,n−1∥ L 2(N−1) N (∂Ω) ∥ui,n∥ L 2(N−1) N−2 (∂Ω) ≤ C̃ ( ∥U∥ L 2(N−1) N−2 (∂Ω) + ∥U∥ L 2(N−1) N−2 (∂Ω) ) ≤ C ′, where C ′ is a constant independent of n. Hence, ∥Un∥(H1(Ω))2 = ∥u1,n∥H1(Ω) + ∥u2,n∥H1(Ω) ≤ C, (2.9) where C is a constant independent of n. Since Un is uniformly bounded in (H1(Ω))2 and (H1(Ω))2 is reflexive, it follows that there exists a subsequence (relabeled) Un which converges weakly to U∗ = (u1,∗, u2,∗) ∈ (H1(Ω))2. Now, let’s show that fi(x, Un) + kui,n converges weakly to fi(x, U∗) + kui,∗ in L 2(N−1) N (∂Ω) for i = 1, 2, where we have adopted the notation fi(x, V ) = fi(x, v1, v2) for V = (v1, v2). From (2.9), we have that the sequence Un is monotone increasing and bounded. Therefore, Un converges pointwise to U∗, that is, U∗(x) = limn→∞ Un(x) and U(x) ≤ U∗(x) ≤ U(x) a.e. x ∈ Ω and a.e. x ∈ ∂Ω. Since fi is continuous with respect to the second and third variables, it follows that lim n→∞ [fi(x, Un(x)) + kui,n(x)] = fi(x, U∗(x)) + kui,∗(x). By the Lebesgue Dominated Convergence Theorem, fi(x, Un(x)) + kui,n(x) converges strongly to fi(x, U∗(x)) + kui,∗ in L 2(N−1) N (∂Ω). Hence, we have fi(x, Un(x)) + kui,n(x) converges weakly to fi(x, U∗(x)) + kui,∗, that is, for all ψ ∈ H1(Ω), lim n→∞ ∫ ∂Ω [fi(x, Un) + kui,n]ψ = ∫ ∂Ω [fi(x, U∗(x)) + kui,∗]ψ. (2.10) 6 S. BANDYOPADHYAY, T. LEWIS, N. MAVINGA EJDE-2025/43 We will show that U∗ = (u1,∗, u2,∗) is a weak solution of (1.1). By the continuity of the trace operator (1.2) and the embedding L 2(N−1) N−2 (∂Ω) into L 2(N−1) N (∂Ω), it follows that ui,∗ ∈ L 2(N−1) N (∂Ω) for i = 1, 2. Using (2.8), we have ∥fi(x, U∗)∥ L 2(N−1) N (∂Ω) = ∥fi(x, U∗) + kui,∗ − kui,∗∥ L 2(N−1) N (∂Ω) ≤ ∥fi(x, U∗) + kui,∗∥ L 2(N−1) N (∂Ω) + ∥kui,∗∥ L 2(N−1) N (∂Ω) ≤ C, where C is a positive constant. Hence, fi(x, U∗) ∈ L 2(N−1) N (∂Ω). Since ui,n converges weakly to ui,∗ in H1(Ω), ui,n converges strongly in L2(∂Ω), and fi(x, Un)+ kui,n converges weakly to fi(x, U∗) + kui,∗ in L 2(N−1) N (∂Ω). By taking the limit in (2.6) as n goes to ∞ and using (2.10), we obtain∫ Ω (∇ui,∗∇ψ + ui,∗ψ) + ∫ ∂Ω kui,∗ψ = lim n→∞ (∫ Ω ( ∇ui,n∇ψ + ui,nψ ) + ∫ ∂Ω kunψ ) = lim n→∞ (∫ ∂Ω (fi(x, Un−1) + kui,n−1)ψ ) = ∫ ∂Ω (fi(x, U∗) + kui,∗)ψ for any ψ ∈ H1(Ω). Hence,∫ Ω (∇ui,∗∇ψ + ui,∗ψ) = ∫ ∂Ω fi(x, U∗)ψ for all ψ ∈ H1(Ω). Thus U∗ is a weak solution to (1.1). Finally, let us show that U∗ is the minimal weak solution to (1.1) in the interval [U,U ]. Let V be a weak solution to (1.1) such that U ≤ V ≤ U . Then V is a weak supersolution of (1.1). Repeating the above iteration procedure with U0 = U , we obtain U ≤ U∗ ≤ V . Thus U∗ is a weak minimal solution. In a similar way, we can construct the maximal weak solution U∗ to (1.1) for which limn→∞Wn = U∗ with W0 = U . This completes the proof of Theorem 1.3. 3. Proof of Theorem 1.4 The proof is based on an application of Zorn’s lemma where we construct a nonempty set of subsolutions and show that the set has a maximal element which will turn out to be a solution to (1.1). We then use a version of Kato’s Inequality up to the boundary to prove that the maximal element of the set of subsolutions is in fact a maximal solution to (1.1). Similarly, we can show the existence of a minimal solution to (1.1) by applying Zorn’s lemma on a set of supersolutions. We proceed to show the existence of a maximal solution. This proof involves several steps described below. Step 1. Existence of a uniformly bounded subsolution of (1.1). Let Ũ = (ũ1, ũ2) ∈ (H(Ω))2 be a subsolution to (1.1) such that U ≤ Ũ ≤ U . We will show the existence of a subsolution W = (w1, w2) of (1.1) such that Ũ ≤ W ≤ U and ∥W∥(H1(Ω))2 ≤ C, where C is a constant depending on ui, ui,Ω,Ki. Consider the equations −∆u1 + u1 = 0 in Ω, ∂u1 ∂η = f1(x, u1, ũ2) on ∂Ω, (3.1) and −∆u2 + u2 = 0 in Ω, ∂u2 ∂η = f2(x, ũ1, u2) on ∂Ω. (3.2) We set f1(x, s1) := f1(x, s1, ũ2(x)). It follows from condition (A4) that whenever ũ1 ≤ s1 ≤ u1, ũ2 ≤ u2 and r > 2(N−1) N , it holds |f1(x, s1)| = |f1(x, s1, ũ2(x))| ≤ K1(x) ∈ Lr(∂Ω). Then, EJDE-2025/43 EXISTENCE OF MAXIMAL AND MINIMAL WEAK SOLUTIONS 7 applying Proposition 1.5 to (3.1), we obtain a solution w1 ∈ H1(Ω) such that ũ1 ≤ w1 ≤ u1. Similarly, setting f2(x, s2) := f2(x, ũ1, s2) and applying Proposition 1.5 to (3.2), we obtain a solution w2 such that ũ2 ≤ w2 ≤ u2. Since fi is quasimonotone non-decreasing, it follows that f1(x,w1, ũ2) ≤ f1(x,w1, w2) and f2(x, ũ1, w2) ≤ f2(x,w1, w2). Hence,∫ Ω (∇wi∇ψ + wiψ) ≤ ∫ ∂Ω fi(x,w1, w2)ψ for all ψ ∈ H1(Ω). Using (A4) and ui ≤ ũi ≤ wi ≤ ui, we have that fi(x,w1, w2) ∈ Lr(∂Ω) with r > 2(N−1) N . Hence, fi(x,w1, w2) ∈ L 2(N−1) N (∂Ω) for i = 1, 2. So, W = (w1, w2) is a subsolution of (1.1). Furthermore, ∥W∥(H1(Ω))2 ≤M(ui, ui,Ω,Ki), where M(ui, ui,Ω,Ki) is a constant. Indeed, by using Hölder’s inequality, (A4) and the continuity of the trace operator, we obtain ∥w1∥2H1(Ω) = ∫ Ω |∇w1|2 + w2 1 = ∫ ∂Ω f1(x,w1, ũ2)w1 ≤ ∫ ∂Ω K1(x)w1 ≤ ∫ ∂Ω K1(x)u1 ≤ ∥K1∥Lr(∂Ω)∥u1∥Lr′ (∂Ω) ≤ C1, where r′ is the conjugate of r and r′ < 2(N−1) N−2 . Similarly ∥w2∥2H1(Ω) ≤ C2. Hence, ∥W∥(H1(Ω))2 ≤ M , where M is constant independent of W . Step 2. Zorn’s Lemma. Consider the set A consisting of (w1, w2) ∈ (H1(Ω))2 such that there exists a subsolution (ũ1, ũ2) of (1.1) satisfying (u1, u2) ≤ (ũ1, ũ2) ≤ (w1, w2) ≤ (u1, u2), where w1 and w2 are solutions of (3.1) and (3.2), respectively, for the pair ũ1, ũ2. We will first check the hypothesis of Zorn’s Lemma, and then we derive the existence of a maximal element of A. From Step 1, we observe that A ̸= ∅. Let Y = {Wn = (w1,n, w2,n)}n≥1 be a chain in (A,≤) where the countability of the indexing set is guaranteed by the separability of the product space (H1(Ω))2. By reordering we have without loss of generality Y = {Wn}n≥1 is an (componentwise) increasing sequence in A. Now, let us show that Y has an upper bound in A. Since Wn = (w1,n, w2,n) belongs to A for every n, there exists a subsolution (ũ1,n, ũ2,n) of (1.1) such that (u1, u2) ≤ (ũ1,n, ũ2,n) ≤ (w1,n, w2,n) ≤ (u1, u2) and w1,n and w2,n are solutions of (3.1) and (3.2), respectively, for the pair ũ1,n, ũ2,n. From Step 1, it follows that ∥Wn∥(H1(Ω))2 ≤ M , where M(ui, ui,Ω,Ki) is a constant independent of n. By the reflexivity of (H1(Ω))2, there is a subsequence (relabeled) Wn which converges weakly to W∗ = (w1,∗, w2,∗). Since the sequence {w1,n} is monotonically increasing and bounded above, {w1,n} converges pointwise to w1,∗. Sim- ilarly, {w2,n} converges pointwise to w2,∗. By the continuity of fi(x, ·, ·), it follows that fi(x,w1,n(x), w2,n(x)) → fi(x,w1,∗(x), w2,∗(x)); i = 1, 2 as n → ∞. Using (A4), |f1(x,w1,n, w2,n)| ≤ K1(x), where K1(x) ∈ Lr(∂Ω) with r > 2(N−1) N−2 > 1. Therefore, by the dominated convergence theorem, we have ∥f1(x,w1,n(x), w2,n(x))− f1(x,w1,∗(x), w2,∗(x))∥Lr(∂Ω) → 0 as n→ ∞. (3.3) Using Holder’s inequality and (3.3), we have that for any test function ψ ∈ H1(Ω),∣∣ ∫ ∂Ω f1(x,w1,n(x), w2,n(x))ψ − ∫ ∂Ω f1(x,w1,∗(x), w2,∗(x))ψ ∣∣ ≤ ∥f1(x,w1,n(x), w2,n(x))− f1(x,w1,∗(x), w2,∗(x))∥Lr(∂Ω).∥ψ∥Lr′ (∂Ω) → 0 as n→ ∞. Therefore, lim n→∞ ∫ ∂Ω f1(x,w1,n(x), w2,n(x))ψ = ∫ ∂Ω f1(x,w1,∗(x), w2,∗(x))ψ (3.4) 8 S. BANDYOPADHYAY, T. LEWIS, N. MAVINGA EJDE-2025/43 for all ψ ∈ H1(Ω). Utilizing (3.4) and the quasimonotonicity condition of f1, we have∫ Ω (∇w1,∗∇ψ + w1,∗ψ) = lim n→∞ ∫ Ω (∇w1,n∇ψ + w1,nψ) = lim n→∞ ∫ ∂Ω f1(x,w1,n(x), ũ2,n(x))ψ ≤ lim n→∞ ∫ ∂Ω f1(x,w1,n(x), w2,n(x))ψ = ∫ ∂Ω f1(x,w1,∗(x), w2,∗(x))ψ and, similarly, ∫ Ω (∇w2,∗∇ψ + w2,∗ψ) ≤ ∫ ∂Ω f2(x,w1,∗(x), w2,∗(x))ψ for all 0 ≤ ψ ∈ H1(Ω). Hence, W∗ = (w1,∗, w2,∗) is a subsolution of (1.1). Taking Ũ = W∗ = (w1,∗, w2,∗) and applying Step 1, there exists a subsolution V = (v1, v2) of (1.1) such that (u1, u2) ≤ (w1,∗, w2,∗) ≤ (v1, v2) ≤ (u1, u2) with v1 and v2 solutions of (3.1), (3.2), respectively, for the pair (w1,∗, w2,∗). This implies V ∈ A and is an upper bound of Y . By Zorn’s Lemma (see e.g. [6]), A has a maximal element Z = (z1, z2) ∈ A. We claim that Z = (z1, z2) is a subsolution of (1.1). Indeed, Z = (z1, z2) ∈ A implies that there exists a subsolution to (1.1), (z̃1, z̃2), such that (u1, u2) ≤ (z̃1, z̃2) ≤ (z1, z2) ≤ (u1, u2) where z1 and z2 are solutions of (3.1) and (3.2), respectively, for the pair z̃1, z̃2. By the quasimonotonicity of f1 and f2, it follows that∫ Ω (∇z1∇ψ + z1ψ) = ∫ ∂Ω f1(x, z1, z̃2)ψ ≤ ∫ ∂Ω f1(x, z1, z2)ψ and ∫ Ω (∇z2∇ψ + z2ψ) = ∫ ∂Ω f2(x, z̃1, z2)ψ ≤ ∫ ∂Ω f2(x, z1, z2)ψ for all 0 ≤ ψ ∈ H1(Ω). Therefore, Z = (z1, z2) is a subsolution of (1.1). Step 3. Z = (z1, z2) is a solution of (1.1). Taking Ũ = Z and applying Step 1, we have that there exists Z∗ = (z∗1 , z ∗ 2) that is a subsolution of (1.1) with zi ≤ z∗i ≤ ui, i = 1, 2, and z∗1 and z∗2 solutions of (3.1) and (3.2), respectively, for the pair (z1, z2). From the definition of the set A, we have that Z∗ = (z∗1 , z ∗ 2) ∈ A. Using the fact that Z is a maximal element of A, we obtain that Z∗ ≤ Z. Hence, Z = Z∗, and so, for every ψ ∈ H1(Ω),∫ Ω (∇z1∇ψ + z1ψ) = ∫ Ω (∇z∗1∇ψ + z∗1ψ) ∫ ∂Ω f1(x, z1, z2)ψ and ∫ Ω (∇z2∇ψ + z2ψ) = ∫ Ω (∇z∗2∇ψ + z∗2ψ) = ∫ ∂Ω f2(x, z1, z2)ψ. Therefore, (z1, z2) is a solution of (1.1) and (u1, u2) ≤ (z1, z2) ≤ (u1, u2). Step 4. Z = (z1, z2) is a maximal solution of (1.1). Let U = (u1, u2) be any solution of (1.1) such that (u1, u2) ≤ (u1, u2) ≤ (u1, u2). Therefore, by Proposition 5.4 in the Appendix (see Section 5), V = (v1, v2) is a subsolution of (1.1) with v1 = max{u1, z1} and v2 = max{u2, z2}. Then it follows from Step 1 that there exists Ẑ = (ẑ1, ẑ2) ∈ A such that U ≤ V ≤ Ẑ ≤ U , which implies V ∈ A. As Z is a maximal element of A, hence, V ≤ Z. On the other hand, since V = (max{u1, z1},max{u2, z2}), we have, Z ≤ V ≤ Ẑ. Consequently, Z = V implying U ≤ Z, and it follows that Z = (z1, z2) is the unique maximal solution of (1.1). By a similar approach, we can show the existence of a minimal solution of (1.1). EJDE-2025/43 EXISTENCE OF MAXIMAL AND MINIMAL WEAK SOLUTIONS 9 4. Finite difference approximations In this section, we use a finite difference (FD) method for approximating solutions of the problem −∆ui + ui = 0 in Ω; ∂ui ∂η = λfi(u1, u2), on ∂Ω, i = 1, 2, (4.1) where λ is parameter and fi : [0,∞) × [0,∞) → R, in addition to (A1), satisfies the following conditions: (i) f1(s1, s2) is locally Lipschitz continuous in s1 and f2(s1, s2) is locally Lipschitz continuous in s2. (ii) fi(0, 0) ≥ 0. (iii) fi’s are sublinear, i.e. lim ∥(s1,s2)∥1→∞ fi(s1, s2) ∥(s1, s2)∥1 = 0, where ∥(s1, s2)∥1 = |s1|+ |s2|. We prove the existence of nonnegative solutions for the discrete problem generated by the FD method (see Section 4.1) in between an ordered pair of discrete sub and supersolutions which turn out to be uniformly bounded independent of the discretization parameter h. This result and the corresponding sub and supersolution technique is a discrete analogue of Theorem 1.3. In fact, we find exact sub and supersolutions for the discrete problem. We formulate a monotone iteration to find the maximal nonnegative solution bounded above by the supersolution. Several bifurcation diagrams generated using MATLAB are provided. 4.1. Formulation. The main idea in our FD formulation is to approximate all differential oper- ators by discrete operators using difference quotients. We first discretize the domain and, at each of the grid points, we approximate the value of the solution by solving the algebraic system of equations that results from replacing the differential operators with discrete difference operators. Assume the domain Ω is an N -rectangle, where N ≥ 1 is the dimension. In other words, Ω = (a1, b1)× (a2, b2)× . . . (aN , bN ). LetMi ≥ 2 be a positive integer and hi = bi−ai Mi−1 for i = 1, 2, . . . N . Define h = (h1, h2, . . . hN ) ∈ RN , M = ∏N i=1(Mi), and NN M = {α = (α1, α2, . . . αN ) | 1 ≤ αi ≤ Mi, i = 1, 2, . . . N}. Next we partition Ω into ∏N i=1(Mi − 1) sub-N rectangles with grid points xα = (a1 + (α1 − 1)h1, a2 + (α2 − 1)h2, . . . , aN + (αN − 1)hN ) for each multi-index α ∈ NN M . We call Th = {xα}α∈NN M a grid for Ω. Let {ei}Ni=1 denote the canonical basis vectors for RN . We define the discrete operators for approximating first order partial derivatives ∂ ∂xi u(x) by δ+xi,hi u(x) := u(x+ hiei)− u(x) hi , δ−xi,hi u(x) := u(x)− u(x− hiei) hi , δxi,hi u(x) := 1 2 δ+xi,hi u(x) + 1 2 δ−xi,hi u(x) = u(x+ hiei)− u(x− hiei) 2hi (4.2) for the function u : RN → R and δ+xi,hi uh(xα) := uh(xα+ei)− uh(xα) hi , δ−xi,hi uh(xα) := uh(xα)− uh(xα−ei) hi , δxi,hi uh(xα) := 1 2 δ+xi,hi uh(xα) + 1 2 δ−xi,hi uh(xα) = uh(xα+ei)− uh(xα−ei) 2hi (4.3) for all xα ∈ Th ∩ Ω for the grid function uh : Th → R. Note that the discrete operators δ±xi,hi are first-order accurate whereas δxi,hi is second-order accurate. We also define the corresponding 10 S. BANDYOPADHYAY, T. LEWIS, N. MAVINGA EJDE-2025/43 discrete gradient operators [∇± h ]i := δ±xi,hi , [∇h]i := δxi,hi . Let ∂̃Ω ⊂ ∂Ω be such that ∂̃Ω := ∂Ω \ {the points where ∂Ω is not smooth}. For x ∈ Th ∩ ∂̃Ω, we define the discrete outer normal derivative using the discrete gradient operator ∇∗ h by [∇∗ hu(x)]i :=  δ+xi,hi u(x) if ηi(x) < 0, δ−xi,hi u(x) if ηi(x) > 0, δxi,hi u(x) if ηi(x) = 0 to ensure that ∇∗ hu(x) · η does not require points outside of the domain Ω. Note that the discrete outward normal derivative approximation is only first order accurate. Next, we define the second order central difference operators for approximating second order nonmixed partial derivatives ∂2 ∂x2 i u(x) by δ2xi,hi u(x) := δ±xi,hi (δ∓xi,hi (u(x))) = u(x+ hiei)− 2u(x) + u(x− hiei) h2i (4.4) for the function u : RN → R and δ2xi,hi uh(xα) := uh(xα+ei)− 2uh(xα) + uh(xα−ei) h2i (4.5) for all xα ∈ Th ∩Ω for the grid function uh : Th → R. Finally, we define the second order discrete Laplacian operator ∆h by ∆h := N∑ i=1 δ2xi,hi . In this section we use the following discrete problem to approximate the solutions to (4.1). Here the grid functions ui,h are an approximation for ui over the grid Th for i = 1, 2: −∆hu1,h + u1,h = 0 in Th ∩ Ω, −∆hu2,h + u2,h = 0 in Th ∩ Ω, ∇∗ hu1,h · η − λf1(u1,h, u2,h) = 0 on Th ∩ ∂̃Ω, ∇∗ hu2,h · η − λf2(u1,h, u2,h) = 0 on Th ∩ ∂̃Ω. (4.6) Remark 4.1. Note that we are eliminating the set of points (a set with measure 0) where the outward normal derivative is not defined. Once Uh is defined over Th∩(Ω∪ ∂̃Ω), it can be extended to Th in post-processing. 4.2. Existence and stability. We use sub and supersolution theory in the discrete setting to prove existence and stability results for solutions to (4.6). First, we define discrete sub and supersolutions of (4.6). We say Uh = (u1,h, u2,h) is a subsolution of (4.6) if it satisfies the following conditions: −∆hu1,h + u1,h ≤ 0 in Th ∩ Ω, −∆hu2,h + u2,h ≤ 0 in Th ∩ Ω, ∇∗ hu1,h · η − λf1(u1,h, u2,h) ≤ 0 on Th ∩ ∂̃Ω, ∇∗ hu2,h · η − λf2(u1,h, u2,h) ≤ 0 on Th ∩ ∂̃Ω. (4.7) We can define the supersolution Uh = (u1,h, u2,h) by reversing the inequalities in (4.7). We focus on proving the existence and stability of nonnegative solutions to (4.6). Step 1. Constructing discrete sub and supersolutions of (4.6). It is clear that φ := (0, 0) is a subsolution of (4.6) since fi(0) ≥ 0 for i = 1, 2. To construct a supersolution of (4.6), let us define EJDE-2025/43 EXISTENCE OF MAXIMAL AND MINIMAL WEAK SOLUTIONS 11 mi := ai+bi 2 for i = 1, 2, . . . , N to be the midpoint of the domain Ω along the xi direction, and, for some c≫ 1 specified later, define the quadratic function ϕ : Ω → R as follows: ϕ(x) := c N∑ i=1 [(xi −mi) 2 + 4]. (4.8) Clearly, ϕ(x) ≥ 4cN on Ω. Observe that ∂2 ∂x2 i ϕ(x) = 2c. Since the operators δ2xi,hi are exact for quadratic functions, we have −∆hϕ(x) + ϕ(x) = −∆ϕ(x) + ϕ(x) = −2cN + c N∑ i=1 (xi −mi) 2 + 4cN = 2cN + c N∑ i=1 (xi −mi) 2 > 0 (4.9) for all x ∈ Th ∩Ω. Choose x ∈ Th ∩ ∂̃Ω, and define Hi = mi − ai = bi−ai 2 . Suppose xi = ai. Then η = −ei, and by the convexity of ϕ, we have ∇∗ hϕ(x) · η = −δ+xi,hi ϕ(x) ≥ −δ+xi,Hi ϕ(x) = −ϕ(x+Hiei)− ϕ(x) Hi = −c(ai +Hi −mi) 2 + 4 + c(ai −mi) 2 − 4 Hi = c bi − ai 2 . Similarly, if xi = bi, then, η = ei, and by the convexity of ϕ(x) it holds ∇∗ hi ϕ(x) · η = δ−xi,hi ϕ(x) ≥ δ−xi,Hi ϕ(x) = ϕ(x)− ϕ(x−Hiei) Hi = c(bi −mi) 2 + 4− c(bi −Hi −mi) 2 − 4 Hi = c bi − ai 2 . Observe that ∥(ϕ(x), ϕ(x))∥1 ≥ 8cN → ∞ as c→ ∞. Furthermore, ∥(ϕ(x), ϕ(x))∥1 ≤ cM, where M = 2( ∑N i=1 4H 2 i + 4N). Hence, by the sublinearity of fi, 0 ≤ fi(ϕ(x), ϕ(x)) ∥(ϕ(x), ϕ(x))∥1 ≤ fi(cM, cM) 8cN = M 8N fi(cM, cM) cM → 0 as c→ ∞ since M 8N is a constant independent of c. Therefore, there exists c≫ 1 such that ∇∗ hϕ(x) · η (∥(ϕ(x), ϕ(x))∥1) ≥ c(bi − ai) 2(∥(ϕ(x), ϕ(x))∥1) ≥ bi − ai ( ∑N k=1 4H 2 k + 4N) > λfj(ϕ(x), ϕ(x)) ∥(ϕ(x), ϕ(x))∥1 (4.10) for j = 1, 2. Thus, ∇∗ hϕ(x) ·η ≥ λfj(ϕ(x), ϕ(x)) for all x ∈ ∂̃Ω and for j = 1, 2. Finally, combining (4.9) and (4.10), we conclude that (ϕ(x), ϕ(x)) is a supersolution of (4.6). 12 S. BANDYOPADHYAY, T. LEWIS, N. MAVINGA EJDE-2025/43 Step 2. Forming a monotone iteration. Let us assume the reaction terms f1, f2 are Lipschitz continuous and let U (0) h be a discrete supersolution of (4.6). Consider the fixed point iteration U (n+1) h = MKU (n) h (4.11) for all n ≥ 0, where K is the maximum of the Lipschitz constants for f1, f2 in (4.1) and MK is defined such that −∆hu (n+1) 1,h + u (n+1) 1,h + λKu (n+1) 1,h = λKu (n) 1,h in Th ∩ Ω, −∆hu (n+1) 2,h + u (n+1) 2,h + λKu (n+1) 2,h = λKu (n) 2,h in Th ∩ Ω, ∇∗ hu (n+1) 1,h · η + λKu (n+1) 1,h = λf1(u (n) 1,h, u (n) 2,h) + λKu (n) 1,h on Th ∩ ∂̃Ω, ∇∗ hu (n+1) 2,h · η + λKu (n+1) 2,h = λf2(u (n) 1,h, u (n) 2,h) + λKu (n) 2,h on Th ∩ ∂̃Ω. (4.12) Before we proceed with the following theorem, let us write the above mapping on grid functions as an equivalent transformation for vectors. Let J0 = |Th ∩ (Ω ∪ (∂̃Ω))| and U ∈ R2J0 denote the vectorization of the grid function Uh. Notationally, the I subscript will correspond to grid function values in Th ∩ Ω and the B subscript will correspond to grid function values in Th ∩ ∂̃Ω. Then, (4.12) is equivalent to MU(n+1) = λF(U(n)), (4.13) where M =  LI 0 LB 0 0 LI 0 LB BI 0 BB 0 0 BI 0 BB  , U =  u1,I u2,I u1,B u2,B  , F =  F1,I F2,I F1,B F2,B  for LI and LB matrices corresponding to −∆h+(1+λK)I; BI and BB matrices corresponding to ∇∗ h · η + λKI; Fi,I corresponding to λK for i = 1, 2; and Fi,B corresponding to λfi(u1,u2) + λK for i = 1, 2. Clearly, M is diagonally dominant since it is positive on the diagonal, non-positive for all off-diagonal terms, and the row sum is positive. Hence, M is a Z-matrix, and by Gershgorin’s Circle Theorem, M is non-singular since the real part of all of its eigenvalues are always positive. Therefore, it follows that M is a monotone matrix. Remark 4.2. Notice that the iteration (4.12) is well defined and fi(s1, s2) +Ksi is increasing in s1 and s2 for i = 1, 2. Theorem 4.3. Let Uh = (u1,h, u2,h) be a nonnegative subsolution of (4.6) and U (0) h = (u (0) 1,h, u (0) 2,h) be a supersolution of (4.6) such that U (0) h ≥ Uh. Then, U (1) h = (u (1) 1,h, u (1) 2,h) = MKU (0) h is a supersolution of (4.6) with Uh ≤ U (1) h ≤ U (0) h . Proof. Observe that, by the definition of MK described in (4.12), for i = 1, 2 in Th ∩ Ω, we have −∆hu (1) i,h + u (1) i,h + λKu (1) i,h = λKu (0) i,h ≥ λKui,h ≥ −∆hui,h + ui,h + λKui,h. (4.14) Also, on Th ∪ ∂̃Ω, we have ∇∗ hu (1) 1,h · η + λKu (1) 1,h = λf1(u (0) 1,h, u (0) 2,h) + λKu (0) 1,h ≥ λf1(u (0) 1,h, u2,h) + λKu (0) 1,h ≥ λf1(u1,h, u2,h) + λKu1,h ≥ ∇∗ hu1,h · η + λKu1,h by the quasimonotonicity of fi and Remark 4.2. In a similar fashion there holds ∇∗ hu (1) 2,h · η + λKu (1) 2,h = λf1(u (0) 1,h, u (0) 2,h) + λKu (0) 2,h ≥ λf1(u1,h, u (0) 2,h) + λKu (0) 2,h ≥ λf1(u1,h, u2,h) + λKu2,h ≥ ∇∗ hu2,h · η + λKu2,h. EJDE-2025/43 EXISTENCE OF MAXIMAL AND MINIMAL WEAK SOLUTIONS 13 Hence MU(1) ≥MU, and it follows that U(1) ≥ U. Thus U (1) h ≥ Uh. Next, for i = 1, 2, −∆hu (1) i,h + u (1) i,h + λKu (1) i,h = λKu (0) i,h = λKu (0) i,h + 0 ≤ −∆hu (0) i,h + u (0) i,h + λKu (0) i,h in Th ∩ Ω since U (0) h is a nonnegative supersolution of (4.6). Furthermore, ∇∗ hu (1) i,h · η + λKu (1) i,h = λf1(u (0) 1,h, u (0) 2,h) + λKu (0) i,h ≤ ∇∗ hu (0) i,h · η + λKu (0) i,h on Th ∩ ∂̃Ω since U (0) h is a nonnegative supersolution of (4.6). Hence, MU(1) ≤ MU(0) which implies U(1) ≤ U(0). Thus, U (1) h ≤ U (0) h . Since, U (1) h ≤ U (0) h , (4.14) implies −∆hu (1) i,h + u (1) i,h = λKu (0) i,h − λKu (1) i,h ≥ 0 in Th ∩ Ω for i = 1, 2. Also, from the definition of MK combined with the quasimonotonicity of fi and Remark 4.2, we obtain ∇∗ hu (1) 1,h · η + λKu (1) 1,h = λf1(u (0) 1,h, u (0) 2,h) + λKu (0) 1,h ≥ λf1(u (1) 1,h, u (0) 2,h) + λKu (1) 1,h ≥ λf1(u (1) 1,h, u (1) 2,h) + λKu (1) 1,h on Th ∩ ∂̃Ω. Thus ∇∗ hu (1) 1,h · η ≥ f1(u (1) 1,h, u (1) 2,h) on Th ∩ ∂̃Ω. Similarly, ∇∗ hu (1) 2,h · η + λKu (1) 2,h = λf2(u (0) 1,h, u (0) 2,h) + λKu (0) 2,h ≥ λf2(u (0) 1,h, u (1) 2,h) + λKu (1) 2,h ≥ λf1(u (1) 1,h, u (1) 2,h) + λKu (1) 2,h on Th ∩ ∂̃Ω. Thus ∇∗ hu (1) 2,h · η ≥ λf2(u (1) 1,h, u (1) 2,h) on Th ∩ ∂̃Ω. Therefore, U (1) h is a supersolution of (4.6). □ Step 3. Existence of a fixed point. Theorem 4.4. Let 0 ≤ Uh be a subsolution to (4.6) and U (0) h ≥ Uh be a supersolution of (4.6). Then the sequence U (n) h defined by (4.12) converges to a solution of (4.6). Furthermore, if ∥U (0) h ∥l∞([Th]2) is bounded independent of h, then the solution is ℓ∞ norm stable. Proof. Observe that, by Theorem 4.3, we have 0 ≤ Uh ≤ U (n+1) h ≤ U (n) h ≤ U (n−1) h ≤ · · · ≤ U (1) h ≤ U (0) h for all n ≥ 1. Thus, the sequence {U (n) h }∞n=0 is convergent since it is monotone and bounded. Let Vh : [T̃h]2 → R2 such that U (n) h → Vh in l∞([T̃h]2) for T̃h = Th ∩ (Ω ∪ ∂̃Ω). Clearly, Vh is a fixed point of (4.11). Thus Vh is a solution of (4.6) with 0 ≤ Uh ≤ Vh ≤ U (n+1) h ≤ U (n) h ≤ U (n−1) h ≤ · · · ≤ U (1) h ≤ U (0) h from which the stability result follows. □ 14 S. BANDYOPADHYAY, T. LEWIS, N. MAVINGA EJDE-2025/43 Remark 4.5. Let ϕ be defined by (4.8) so that it is bounded by a constant depending only on Ω, λ, and f for f satisfying the hypotheses for problem (4.1). Then ϕ is bounded independent of h, and it follows that the grid functions Uh defined by Uh(x) := (0, 0) for all x ∈ Th and U (0) h defined by U (0) h (x) := (ϕ(x), ϕ(x)) for all x ∈ Th for c sufficiently large satisfy the assumptions in Theorem 4.4 to ensure the finite difference approximation exists and is ℓ∞ norm stable. 4.3. Example 1. In this example, we consider the nonlinearities f1(u1, u2) = 100 √ u2 − 2 √ u1, f2(u1, u2) = 4 √ u1 + √ u2. (4.15) Clearly the fi’s are Caratheódory, quasimonotone non-decreasing and sublinear functions satisfying the assumption (A2) and (A3) in Theorem 1.3. Moreover, fi(0, 0) = 0, ∂f1 ∂u2 (0, 0) = ∞ and ∂f2 ∂u1 (0, 0) = ∞. Solutions can be found in Figure 1 for λ = 0.5 and λ = 3, and computed bifurcation diagrams can be found in Figure 2. When choosing the supersolution U (0) h according to Remark 4.5, the solution for λ = 0.5 was found using c = 16, 384 and the solution for λ = 3 was found using c = 524, 288. Figure 1. Graphs of U1 and U2 for λ = 0.5 and λ = 3 Figure 2. Bifurcation diagrams of (4.1) for λ ∈ (0, 10) when f1(u1, u2) = 100 √ u2 − 2 √ u1, f2(u1, u2) = 4 √ u1 + √ u2. 4.4. Example 2. In this example, we consider the nonlinearities f1(u1, u2) = arctan(u2), f2(u1, u2) = 3 arctan(u1). (4.16) EJDE-2025/43 EXISTENCE OF MAXIMAL AND MINIMAL WEAK SOLUTIONS 15 Clearly the fi’s are Caratheodory, quasimonotone, non-decreasing, and sublinear functions satis- fying the assumption (A2) and (A3) in Theorem 1.3. Moreover, fi(0, 0) = 0, ∂f1 ∂u2 (0, 0) > 0 and ∂f2 ∂u1 (0, 0) > 0. Solutions can be found in Figure 3 for λ = 0.5 and λ = 3, and computed bifurcation diagrams can be found in Figure 4. We see nonexistence of a positive solution for λ small. When choosing the supersolution U (0) h according to Remark 4.5, the solution for λ = 0.5 was found using c = 4 and the solution for λ = 3 was found using c = 16. Figure 3. Graphs of U1 and U2 for λ = 0.5 and λ = 3. Bifurcation diagram when λ ∈ (0, 1). Bifurcation diagram when λ ∈ (0, 10). Figure 4. Bifurcation diagrams of (4.1) for different ranges of λ when f1(u1, u2) = arctan(u2), f2(u1, u2) = 3 arctan(u1). 4.5. Example 3. In this example, we consider the nonlinearities f1(u1, u2) = e 10u1 1+u2 1 − 1 + 5 3 √ u22 + 1, f2(u1, u2) = e u2 1+u2 2 − 1. Clearly the fi’s are Caratheodory, quasimonotone, non-decreasing, and sublinear functions sat- isfying the assumption (A2) and (A3) in Theorem 1.3. Moreover, fi(0, 0) = 0, ∂f1 ∂u2 (0, 0) = ∞ and ∂f2 ∂u1 (0, 0) = 0. Solutions can be found in Figure 5 for λ = 0.5 and λ = 3, and computed bifurcation diagrams can be found in Figure 6. Also in Figure 6, notice that there exists a range of λ for which we see non coexistence of positive solutions. Our method and monotone solver found the maximal solution based on our supersolution. Observe that, in Subfigure 6(A), the vertical lines are jump discontinuities corresponding to finding maximal solutions. We used the method of 16 S. BANDYOPADHYAY, T. LEWIS, N. MAVINGA EJDE-2025/43 continuation to find the exact branches of the bifurcation curves in Subfigure 6(B). When choosing the supersolution U (0) h according to Remark 4.5, the solution for λ = 0.5 was found using c = 512 and the solution for λ = 3 was found using c = 65536. Figure 5. Graphs of U1 and U2 for λ = 0.5 and λ = 3. Bifurcation diagram when λ ∈ (0, 0.5). Bifurcation diagram when λ ∈ (0, 1). Figure 6. Bifurcation diagrams of (4.1) for different ranges of λ when f1(u1, u2) = e 10u1 1+u2 1 − 1 + 5 3 √ u22 + 1, f2(u1, u2) = e u2 1+u2 2 − 1. 4.6. Example 4. In this example, we consider the nonlinearities f1(u1, u2) = e u1 1+u2 1 − 1 + 1 3 3 √ u22 + 1, f2(u1, u2) = arctan(u2) + e u2 1+u2 2 − 1. Clearly the fi’s are Caratheodory, quasimonotone, non-decreasing, and sublinear functions satis- fying the assumption (A2) and (A3) in Theorem 1.3. Moreover, fi(0, 0) = 0, ∂f1 ∂u2 (0, 0) = ∞ and ∂f2 ∂u1 (0, 0) > 0. Solutions can be found in Figure 7 for λ = 0.5 and λ = 3, and computed bifurcation diagrams can be found in Figure 8. In Figure 7, notice that for λ = 0.5, u2 > u1, whereas, for λ = 3, u1 > u2. When choosing the supersolution U (0) h according to Remark 4.5, the solution for λ = 0.5 was found using c = 1 and the solution for λ = 3 was found using c = 32. EJDE-2025/43 EXISTENCE OF MAXIMAL AND MINIMAL WEAK SOLUTIONS 17 Figure 7. Graphs of U1 and U2 for λ = 0.5 and λ = 3. Bifurcation diagram when λ ∈ (0, 0.5). Bifurcation diagram when λ ∈ (0, 10). Figure 8. Bifurcation diagrams of (4.1) for different ranges of λ when f1(u1, u2) = e u1 1+u2 1 − 1 + 1 3 3 √ u22 + 1, f2(u1, u2) = arctan(u2) + e u2 1+u2 2 − 1. 5. Appendix In this section we use Kato’s inequality up to the boundary for single equations to prove that the componentwise maximum of two solutions of (1.1) is a subsolution of (1.1) (see Proposition 5.4) and componentwise minimum of two solutions of (1.1) is a supersolution of (1.1) (see Proposition 5.6), which we used in Step 4 of the proof of Theorem 1.4 (see Section 3). Now we have Kato’s inequality up to the boundary for single equations. Proposition 5.1. Let u1 and u2 be functions in H1(Ω) such that there exist f1 and f2 in Lr(∂Ω), for r ≥ 2(N−1) N , satisfying∫ Ω (∇ui∇ψ + uiψ) ≤ ∫ ∂Ω fiψ for all 0 ≤ ψ ∈ H1(Ω) , (5.1) for i = 1, 2. Then, u := max{u1, u2} satisfies∫ Ω (∇u∇ψ + uψ) ≤ ∫ ∂Ω fψ for all 0 ≤ ψ ∈ H1(Ω) , where f(x) := { f1(x) if u1(x) > u2(x) a.e. x ∈ ∂Ω, f2(x) if u1(x) ≤ u2(x) a.e. x ∈ ∂Ω For a proof of the above proposition, see [4, Theorem 2.4]. 18 S. BANDYOPADHYAY, T. LEWIS, N. MAVINGA EJDE-2025/43 Corollary 5.2. Let u1 and u2 be functions in H1(Ω) such that there exist f1 and f2 in Lr(∂Ω), for r ≥ 2(N−1) N , satisfying∫ Ω (∇ui∇ψ + uiψ) ≥ ∫ ∂Ω fiψ for all 0 ≤ ψ ∈ H1(Ω) , for i = 1, 2. Then, u := min{u1, u2} satisfies∫ Ω (∇u∇ψ + uψ) ≥ ∫ ∂Ω fψ, for all 0 ≤ ψ ∈ H1(Ω) , where f(x) := { f1(x) if u1(x) < u2(x) a.e. x ∈ ∂Ω f2(x) if u1(x) ≥ u2(x) a.e. x ∈ ∂Ω. For a proof of the above corollary, see [4, Corollary 2.5]. Let A be a set consisting of (w1, w2) ∈ (H1(Ω))2 such that there exists a subsolution (ũ1, ũ2) of (1.1) satisfying (u1, u2) ≤ (ũ1, ũ2) ≤ (w1, w2) ≤ (u1, u2), where w1 and w2 are solutions of (3.1) and (3.2), respectively, for the pair (ũ1, ũ2). Lemma 5.3. Let (α1, α2) and (β1, β2) ∈ A be any two subsolutions of (1.1). Then, the pair (max{α1, β1},max{α2, β2}) is a subsolution of (1.1). Proof. Since (α1, α2) and (β1, β2) belong to A, there exist (α̃1, α̃2) and (β̃1, β̃2) such that ui ≤ α̃i ≤ αi ≤ ui, i = 1, 2, and ui ≤ β̃i ≤ βi ≤ ui, i = 1, 2, which satisfy −∆α1 + α1 = 0 in Ω, ∂α1 ∂η = f1(x, α1, α̃2) on ∂Ω, (5.2) −∆α2 + α2 = 0 in Ω, ∂α2 ∂η = f2(x, α̃1, α2) on ∂Ω, (5.3) −∆β1 + β1 = 0 in Ω, ∂β1 ∂η = f1(x, β1, β̃2) on ∂Ω, (5.4) −∆β2 + β2 = 0 in Ω, ∂β2 ∂η = f2(x, β̃1, β2) on ∂Ω. (5.5) We define γ1 := max{α1, β1} and γ2 := max{α2, β2}. By the quasimonotonicity of fi, the following inequalities hold: f1(x, α1, α̃2) ≤ f1(x, α1, γ2), f1(x, β1, β̃2) ≤ f1(x, β1, γ2), f2(x, α̃1, α2) ≤ f2(x, γ1, α2), f2(x, β̃1, β2) ≤ f2(x, γ1, β2). Let g1(x) := { f1(x, α1, γ2) if α1(x) > β1(x) a.e. x ∈ ∂Ω, f1(x, β1, γ2) if α1(x) ≤ β1(x) a.e. x ∈ ∂Ω and g2(x) := { f2(x, γ1, α2) if α2(x) > β2(x) a.e. x ∈ ∂Ω, f2(x, γ1, β2) if α2(x) ≤ β2(x) a.e. x ∈ ∂Ω. EJDE-2025/43 EXISTENCE OF MAXIMAL AND MINIMAL WEAK SOLUTIONS 19 Notice that α1 and β1 are subsolutions of −∆v + v = 0 in Ω, ∂v ∂η = g1(x) on ∂Ω, (5.6) and α2 and β2 are subsolutions of −∆v + v = 0 in Ω, ∂v ∂η = g2(x) on ∂Ω. (5.7) Then, by Proposition 5.1, along with equations (5.6) and (5.7), we have that (γ1, γ2) is a subsolu- tion of (1.1). Thus, the pair (max{α1, β1},max{α2, β2}) is a subsolution of (1.1) which completes the proof. □ Proposition 5.4. Suppose that (u1, u2) and (v1, v2) are two solutions of (1.1) such that (u1, u2) ≤ (u1, u2), (v1, v2) ≤ (u1, u2), where (u1, u2) and (u1, u2) are sub- and supersolutions of (1.1), re- spectively. Then (max{u1, v1},max{u2, v2}) is a subsolution of (1.1). Proof. Observe that any solution (u1, u2) of (1.1) such that (u1, u2) ≤ (u1, u2) ≤ (u1, u2) belongs to A as (u1, u2) is a solution of (3.1) and (3.2) for the pair (u1, u2). The rest of the proof follows from Lemma 5.3. □ Next, we show that if (u1, u2) and (v1, v2) are two solutions of (1.1), then the pair of functions (min{u1, v1},min{u2, v2}) is a supersolution of (1.1). For that purpose, consider the equations −∆u1 + u1 = 0 in Ω, ∂u1 ∂η = f1(x, u1, û2) on ∂Ω, (5.8) and −∆u2 + u2 = 0 in Ω, ∂u2 ∂η = f2(x, û1, u2) on ∂Ω. (5.9) Let B consisting of (w1, w2) ∈ (H1(Ω))2 such that there exists a supersolution (û1, û2) of (1.1) satisfying (u1, u2) ≤ (w1, w2) ≤ (û1, û2) ≤ (u1, u2), (5.10) where w1 and w2 are solutions of (5.8) and (5.9), respectively, for the pair û1, û2. Lemma 5.5. Let (α1, α2) and (β1, β2) ∈ B be any two supersolutions of (1.1). Then the pair (min{α1, β1},min{α2, β2}) is a supersolution of (1.1). Proof. Since (α1, α2) and (β1, β2) belong to B, there exist (α̂1, α̂2) and (β̂1, β̂2) such that ui ≤ αi ≤ α̂i ≤ ui, i = 1, 2 and ui ≤ βi ≤ β̂i ≤ ui, i = 1, 2, which satisfy −∆α1 + α1 = 0 in Ω, ∂α1 ∂η = f1(x, α1, α̂2) on ∂Ω, (5.11) −∆α2 + α2 = 0 in Ω, ∂α2 ∂η = f2(x, α̂1, α2) on ∂Ω, (5.12) −∆β1 + β1 = 0 in Ω ∂β1 ∂η = f1(x, β1, β̂2) on ∂Ω, (5.13) −∆β2 + β2 = 0 in Ω, ∂β2 ∂η = f2(x, β̂1, β2) on ∂Ω. (5.14) 20 S. BANDYOPADHYAY, T. LEWIS, N. MAVINGA EJDE-2025/43 Now, define γ1 := min{α1, β1} and γ2 := min{α2, β2}. Observe that the quasimonotonicity of fi leads to the inequalities f1(x, α1, α̂2) ≥ f1(x, α1, γ2), f1(x, β1, β̂2) ≥ f1(x, β1, γ2) (5.15) f2(x, α̂1, α2) ≥ f2(x, γ1, α2), f2(x, β̂1, β2) ≥ f2(x, γ1, β2). (5.16) Let g1(x) := { f1(x, α1, γ2) if α1(x) < β1(x) a.e. x ∈ ∂Ω, f1(x, β1, γ2) if α1(x) ≥ β1(x) a.e. x ∈ ∂Ω, and g2(x) := { f2(x, γ1, α2) if α2(x) < β2(x) a.e. x ∈ ∂Ω, f2(x, γ1, β2) if α2(x) ≥ β2(x) a.e. x ∈ ∂Ω. Notice that α1 and β1 are supersolutions of −∆v + v = 0 in Ω, ∂v ∂η = g1(x) on ∂Ω, (5.17) and α2 and β2 are supersolutions of the following: −∆v + v = 0 in Ω, ∂v ∂η = g2(x) on ∂Ω. (5.18) Then, by Corollary 5.2, along with equations (5.17) and (5.18), we have that (γ1, γ2) is a super- solution of (1.1). This completes the proof. □ Proposition 5.6. Suppose that (u1, u2) and (v1, v2) are two solutions of (1.1) such that (u1, u2) ≤ (u1, u2), (v1, v2) ≤ (u1, u2), where (u1, u2) and (u1, u2) are sub- and supersolutions of (1.1), respectively. Then (min{u1, v1},min{u2, v2}) is a supersolution of (1.1). Proof. Observe that any solution (u1, u2) of (1.1) such that (u1, u2) ≤ (u1, u2) ≤ (u1, u2) belongs to B. The rest of the proof follows from Lemma 5.5. □ Acknowledgments. S. Bandyopadhyay and T. Lewis were partially supported by the National Science Foundation (NSF) grant DMS-2111059. N. Mavinga was supported by the Simons Foun- dation grant 965180. This material is also based upon work partially supported by the NSF grant 1440140 which helped support the first and third authors’ residence at the Mathematical Sciences Research Institute in Berkeley, California, during the summer of 2022. References [1] Robert A. Adams, John J. F. Fournier; Sobolev spaces, second ed., Pure and Applied Mathematics (Amster- dam), vol. 140, Elsevier/Academic Press, Amsterdam, 2003. MR: 2424078 [2] Herbert Amann; Nonlinear elliptic equations with nonlinear boundary conditions, New developments in differ- ential equations (Proc. 2nd Scheveningen Conf., Scheveningen, 1975), North-Holland Math. Studies, Vol. 21, North-Holland, Amsterdam, 1976, pp. 43–63. MR: 0509487 [3] P. Amster, P. De Nápoli; A quasilinearization method for elliptic problems with a nonlinear boundary condition, Nonlinear Anal. 66 (2007), no. 10, 2255–2263. MR: 2311029 [4] S. Bandyopadhyay, M. Chhetri, B. B. Delgado, N. Mavinga, R. Pardo; Maximal and minimal weak solutions for elliptic problems with nonlinearity on the boundary, Electron. Res. Arch. 30 (2022), no. 6, 2121–2137. MR: 4410308 [5] Robert Stephen Cantrell, Chris Cosner; On the effects of nonlinear boundary conditions in diffusive logistic equations on bounded domains, J. Differential Equations 231 (2006), no. 2, 768–804. MR: 2287906 [6] Siegfried Carl, Vy Khoi Le, Dumitru Motreanu; Nonsmooth variational problems and their inequalities, Springer Monographs in Mathematics, Springer, New York, 2007, Comparison principles and applications. MR: 2267795 EJDE-2025/43 EXISTENCE OF MAXIMAL AND MINIMAL WEAK SOLUTIONS 21 [7] Mabel Cuesta Leon; Existence results for quasilinear problems via ordered sub- and supersolutions, Ann. Fac. Sci. Toulouse Math. (6) 6 (1997), no. 4, 591–608. MR: 1624298 [8] Lawrence C. Evans; Partial differential equations, second ed., Graduate Studies in Mathematics, vol. 19, American Mathematical Society, Providence, RI, 2010. MR: 2597943 [9] Peter V. Gordon, Eunkyung Ko, R. Shivaji; Multiplicity and uniqueness of positive solutions for elliptic equations with nonlinear boundary conditions arising in a theory of thermal explosion, Nonlinear Anal. Real World Appl. 15 (2014), 51–57. MR: 3110553 [10] D. B. Gurung, K. C. Gokul, P. R. Adhikary; Mathematical model of thermal effects of blinking in human eye, Int. J. Biomath. 9 (2016), no. 1, 1650006, 24. MR: 3417712 [11] Friedhelm Inkmann; Existence and multiplicity theorems for semilinear elliptic equations with nonlinear bound- ary conditions, Indiana Univ. Math. J. 31 (1982), no. 2, 213–221. MR: 648172 [12] Alois Kufner, Oldřich John, Svatopluk Fuč́ık; Function spaces, Noordhoff International Publishing, Leyden; Academia, Prague, 1977, Monographs and Textbooks on Mechanics of Solids and Fluids; Mechanics: Analysis. MR: 0482102 [13] A. A. Lacey, J. R. Ockendon, J. Sabina; Multidimensional reaction diffusion equations with nonlinear boundary conditions, SIAM J. Appl. Math. 58 (1998), no. 5, 1622–1647. MR: 1637882 [14] Olga A. Ladyzhenskaya, Nina N. Ural’tseva; Linear and quasilinear elliptic equations, Translated from the Russian by Scripta Technica, Inc. Translation editor: Leon Ehrenpreis, Academic Press, New York-London, 1968. MR: 0244627 [15] T. Lewis, Q. Morris, Y. Zhang; Convergence, stabilty analysis, and solvers for approximating sublinear positone and semipositone boundary value problems using finite difference method, J. Comp. Appl. Math. 404 (2022), 377–427. [16] Ping Liu, Junping Shi; Bifurcation of positive solutions to scalar reaction-diffusion equations with nonlinear boundary condition, J. Differential Equations 264 (2018), no. 1, 425–454. MR: 3712947 [17] Nsoki Mavinga, Rosa Pardo; Bifurcation from infinity for reaction-diffusion equations under nonlinear bound- ary conditions, Proc. Roy. Soc. Edinburgh Sect. A 147 (2017), no. 3, 649–671. MR: 3656708 [18] Enzo Mitidieri, Guido Sweers; Existence of a maximal solution for quasimonotone elliptic systems, Differential Integral Equations 7 (1994), no. 5-6, 1495–1510. MR: 1269667 [19] J. M. Neuberger, N. Sieben, J. W. Swift; Symmetry and automated branch following for a semilinear elliptic pde on a fractal region, SIAM J. Appl. Dyn. Syst. 5 (2006), 476–507. [20] J. M. Neuberger, J.W. Swift; Newton’s method and morse index for semilinear elliptic pdes, Internat. J. Bifur. Chaos Appl. Sci. Engrg. 11 (2001), 801–820. [21] H. Othmer, R. Dillon; A mathematical model for outgrowth and spatial patterning of the vertebrate limb bud, J. Theoret. Biol. 197 (1999), 295–330. MR: 3918372 [22] C. V. Pao; Nonlinear parabolic and elliptic equations, Plenum Press, New York, 1992. MR: 1212084 [23] C. V. Pao; Monotone iterative methods for finite difference system of reaction-diffusion equations, Numer. Math. 46(4) (1985), 571–586. [24] C. V. Pao; Comparison methods and stabiity analysis of reaction diffusion systems, in: Comparison methods and stability theory, CRC Press (1994), 277–292. [25] C. V. Pao; Block monotone iterative methods for numerical solutions of nonlinear elliptic equations, Numer. Math. 72 (1995), 239–261. [26] Kenichiro Umezu; Uniqueness of a positive solution for the Laplace equation with indefinite superlinear bound- ary condition, J. Differential Equations 350 (2023), 124–151. MR: 4528560 Shalmali Bandyopadhyay Department of Mathematics and Statistics, University of Tennessee at Martin, Martin, TN 38238, USA Email address: sbandyo5@utm.edu Thomas Lewis Department of Mathematics and Statistics, The University of North Carolina at Greensboro, Greens- boro, NC 27412, USA Email address: tllewis3@uncg.edu Nsoki Mavinga Department of Mathematics and Statistics, Swarthmore College, Swarthmore, PA 19081, USA Email address: nmaving1@swarthmore.edu 1. Introduction 2. Proof of Theorem ?? 3. Proof of Theorem ?? 4. Finite difference approximations 4.1. Formulation 4.2. Existence and stability 4.3. Example 1 4.4. Example 2 4.5. Example 3 4.6. Example 4 5. Appendix Acknowledgments References