\documentclass[reqno]{amsart} \AtBeginDocument{{\noindent\small {\em Electronic Journal of Differential Equations}, Vol. 2004(2004), No. 93, pp. 1--15.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu \newline ftp ejde.math.txstate.edu (login: ftp)} \thanks{\copyright 2004 Texas State University - San Marcos.} \vspace{9mm}} \begin{document} \title[\hfilneg EJDE-2004/93\hfil Existence of solutions] {Existence of solutions for a nonlinear degenerate elliptic system} \author[N. M. Chuong, T. D. Ke\hfil EJDE-2004/93\hfilneg] {Nguyen Minh Chuong, Tran Dinh Ke } % in alphabetical order \address{Nguyen Minh Chuong \hfill\break Institute of Mathematics, Vietnam Academy of Science and Technology, 18 Hoang Quoc Viet, 10307 Hanoi, Vietnam} \email{nmchuong@math.ac.vn} \address{Tran Dinh Ke \hfill\break Faculty of Mathematics, Hanoi University of Education, Hanoi, Vietnam} \email{ketd@hn.vnn.vn} \date{} \thanks{Submitted February 18, 2004. Published July 27, 2004.} \thanks{Supported by the National Fundamental Research Program, Vietnam Academy of Science \hfill\break\indent and Technology.} \subjclass[2000]{34B27, 35J60, 35B05} \keywords{Subsolution; supersolution; degenerate elliptic equation} \begin{abstract} In this paper, we study the existence of solutions for degenerate elliptic systems. We use the sub-super solution method, and the existence of classical and weak solutions. Some sub-supersolutions are constructed explicitly, when the nonlinearities have critical or supercritical growth. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{theorem}{Theorem}[section] \newtheorem{lemma}[theorem]{Lemma} \newtheorem{proposition}[theorem]{Proposition} \newtheorem{corollary}[theorem]{Corollary} \section{Introduction} Let $\Omega$ be a smooth bounded domain in $\mathbb{R}^{N_1} \times \mathbb{R}^{N_2}$ with boundary $\partial \Omega$, and having $\{0\}\in \Omega$. We consider the Dirichlet problem $$\begin{gathered} \Delta_x u + |x|^{2k} \Delta_y u + f(x,y,u,v) = 0 \quad \text{in } \Omega, \\ \Delta_x v + |x|^{2k} \Delta_y v + g(x,y,u,v) = 0 \quad \text{in } \Omega, \\ u = v = 0 \quad \text{on } \partial \Omega, \end{gathered} \label{e1}$$ where $\Delta_x = \sum_{i=1}^{N_1} \frac {\partial^2} {\partial x^2}$, $\Delta_y = \sum_{i=1}^{N_2} \frac {\partial^2} {\partial y^2}, k\ge 0, f$ and $g$ are real-valued functions defined on $\Omega \times \mathbb{R}^2$, satisfying certain conditions which will be specified in next sections. We assume in this paper that $N_1, N_2 \ge 1$ and $N(k)=N_1 + (k+1)N_2\ge 3$. Let $\nu=(\nu_x, \nu_y)$ be the outward unit normal to $\partial \Omega$. When the degenerating factor is removed (i.e. $k=0$), system \eqref{e1} reduces to a problem with the Laplace operator. Such systems have been the subject for many studies. In almost all of them, the systems are in Hamiltonian and potential forms and are considered by using variational methods (see \cite{a1,b1} and references there in). The operator $G_k = \Delta_x+|x|^{2k}\Delta_y$ is of a Grushin type which was studied in \cite{g1, t2}. In particular, existence results for problem $$\begin{gathered} \Delta_x u + |x|^{2k} \Delta_y u + f(u) = 0 \quad \text{in } \Omega, \\ u = 0 \quad \text{on } \partial \Omega, \\ \end{gathered} \label{e2}$$ are obtained in \cite{t1}. Moreover, the authors proved the Sobolev embedding theorem and set the critical exponent to $\frac{N(k)+2}{N(k)-2}$. In \cite{k1}, we introduced an existence result for the Hamiltonian system with $G_k$ involving nonlinerities that may change signs. In this work, we use sub-super solutions method to get both classical and weak solutions, even when the nonlinearities have critical or supercritical growth. In the next section, we extend the nonexistence results in \cite{c2} to our system when it has the Hamiltonian form. The study is based on our generalized Pohozaev identity. Section 3 shows the existence of classical solutions for \eqref{e1}. Finally, in section 4, we construct the maximal and minimal weak solutions of our problem. \section{Nonexistence results} In this section, we prove nonexistence results when the domain has a special shape described in the following definition. \noindent {\bf Definition.} A domain $\Omega$ is called k-starshaped with respect to $\{0,0\}$ if the inequality $$(x, \nu_x) + (k+1)(y, \nu_y)\ge 0$$ holds almost everywhere on $\partial \Omega$. \smallskip We are now in position to build a generalized version of Pohozaev identity. Firstly, note that the Pohozaev identity for potential system with Laplace operator was introduced in \cite{a1}. Recall also that, the similar identity for scalar case was obtained in \cite{c2}. Precisely, let $u(x,y)\in H^2(\Omega)$ (the usual Sobolev space) be a solution of the problem $$\begin{gathered} \Delta_x u + |x|^{2k} \Delta_y u + f(u) = 0 \quad \text{in } \Omega, \\ u = 0 \quad \text{on } \partial \Omega, \\ \end{gathered} \label{e2.1}$$ then \begin{align*} &N(k) \int_{\Omega} F(u) \,dx\,dy - \frac{N(k) - 2}2 \int_{\Omega} f(u)u \,dx\,dy \\ &=\frac 12 \int_{\partial \Omega} [(x, \nu_x) + (k+1)(y, \nu_y)] (|\nu_x|^2 + |x|^{2k}|\nu_y|^2)\Big( \frac {\partial u}{\partial \nu}\Big)^2 dS, \end{align*} where $F(u) = \int_0^u f(s) ds$. In our paper, Lemma \ref{lm1} makes an extension of this identity for the Hamiltonian system with Grushin operator. Before stating Lemma \ref{lm1}, we state the following condition. \begin{itemize} \item{(S1)} There exists a function $H(x, y, u, v) \in C^1(\Omega \times \mathbb{R}^2)$ satisfying $$\frac {\partial H}{\partial v} = f, \quad \frac {\partial H}{\partial u} = g, \quad H(x,y,0,0)=0 \quad \text{for } (x,y)\in \Omega.$$ \end{itemize} For the conditions in (S1), one can take $f=f(v), g=g(u)$. Thus, $H(u,v) = F(v) + G(u)$, where $$F(v) = \int_0^v f(s) ds\,,\quad G(u) = \int_0^u g(t) dt\,.$$ Denote $$\nabla_x = (\frac {\partial}{\partial x_1},\dots, \frac {\partial}{\partial x_{N_1}}), \quad \nabla_y = (\frac {\partial}{\partial y_1},\dots,\frac {\partial}{\partial y_{N_2}}).$$ \begin{lemma}[Generalized Pohozaev identity] \label{lm1} Let $\Omega$ be a k-starshaped domain with respect to $\{0,0\}$ and let (S1) hold. If $(u,v)\in H^2(\Omega)\times H^2(\Omega)$ is a solution of \eqref{e1} then $(u,v)$ satisfies the equation \begin{align*} &N(k) \int_{\Omega} H(x,y,u,v) \,dx\,dy + \int_{\Omega} [(x,\nabla_x H) + (k+1)(y,\nabla_y H)]dx\,dy \\ &=(N(k)-2) \int_{\Omega} [tf(x,y,u,v)v + (1-t)g(x,y,u,v)u] dx\,dy \\ &\quad + \int_{\partial \Omega} [(x,\nu_x) + (k+1)(y, \nu_y)](|\nu_x|^2 + |x|^{2k}|\nu_y|^2)\frac {\partial u}{\partial \nu}\frac{\partial v}{\partial \nu} dS, \end{align*} for all $t\in \mathbb{R}$. \end{lemma} \begin{proof} For $i=1,\dots ,N_1$, we have \begin{align*} &\int_{\Omega} \frac{\partial}{\partial x_i} \Big(\ x_i H(x,y,u,v)\Big) dy\,dy \\ &= \int_{\Omega} H(x,y,u,v) \,dx\,dy + \int_{\Omega} x_i (g\frac{\partial u}{\partial x_i} + f\frac{\partial v}{\partial x_i} + \frac{\partial H}{\partial x_i}) \,dx\,dy = 0. \end{align*} This implies $$\int_{\Omega} H(x,y,u,v) \,dx\,dy = -\int_{\Omega} x_i [g\frac{\partial u}{\partial x_i} + f\frac{\partial v}{\partial x_i} + \frac {\partial H}{\partial x_i}]\,dx\,dy.$$ Hence $$\int_{\Omega} H(x,y,u,v) \,dx\,dy = -\frac 1{N_1} \int_{\Omega} [(x, \nabla_x u)g + (x, \nabla_x v)f + (x, \nabla_x H)]\,dx\,dy\,. \label{e3}$$ Analogously, for $\beta \in \mathbb{R}$, $$\beta \int_{\Omega} H(x,y,u,v) \,dx\,dy = -\frac {\beta}{N_2} \int_{\Omega} [(y, \nabla_y u)g + (y, \nabla_y v)f + (y, \nabla_y H)]\,dx\,dy. \label{e4}$$ Equalities \eqref{e3} and \eqref{e4} yield \begin{aligned} (1+\beta) \int_{\Omega} H(x,y,u,v) \,dx\,dy &= \int_{\Omega} \Big[\frac{(x, \nabla_x u)}{N_1} + \frac{\beta (y, \nabla_y u)}{N_2}\Big]G_k v \,dx\,dy \\ &\quad + \int_{\Omega} \Big[\frac{(x, \nabla_x v)}{N_1} + \frac{\beta (y, \nabla_y v)}{N_2}\Big]G_k u \,dx\,dy \\ &\quad - \int_{\Omega} \Big[\frac{(x, \nabla_x H)}{N_1} + \frac{\beta (y, \nabla_y H)}{N_2}\Big]\,dx\,dy\,. \end{aligned}\label{e5} We make some computations for the following integrals \begin{gather*} I_1 = \frac 1{N_1} \int_{\Omega} [(x, \nabla_x u)\Delta_x v + (x, \nabla_x v) \Delta_x u] \,dx\,dy, \\ I_2 = \frac {\beta}{N_2} \int_{\Omega} [(y, \nabla_y u)\Delta_x v + (y, \nabla_y v) \Delta_x u] \,dx\,dy, \\ I_3 = \frac 1{N_1} \int_{\Omega} [(x, \nabla_x u)\Delta_y v + (x, \nabla_x v) \Delta_y u]|x|^{2k} \,dx\,dy, \\ I_4 = \frac {\beta}{N_2} \int_{\Omega} [(y, \nabla_y u)\Delta_y v + (y, \nabla_y v) \Delta_y u]|x|^{2k} \,dx\,dy. \end{gather*} We have the generalized Rellich identity (for detail computations, we refer the reader to \cite{c2}) \begin{aligned} I_1 &= \frac 1{N_1} \int_{\Omega} [(x, \nabla_x u)\Delta_x v + (x, \nabla_x v)\Delta_x u]\,dx\,dy \\ &= \frac {N_1 -2}{N_1} \int_{\Omega} \nabla_x u\nabla_x v \,dx\,dy + \frac 1{N_1} \int_{\partial \Omega} |\nu_x|^2 (x, \nu_x) \frac {\partial u}{\partial \nu}\frac{\partial v}{\partial \nu}dS. \end{aligned} \label{e10} In the same way \begin{aligned} I_2 &= \frac {\beta}{N_2} \int_{\Omega} [(y, \nabla_y u)\Delta_x v + (y, \nabla_y v)\Delta_x u]\,dx\,dy \\ &= \beta \int_{\Omega} \nabla_x u\nabla_x v \,dx\,dy + \frac {\beta}{N_2} \int_{\partial \Omega} |\nu_x|^2 (y, \nu_y) \frac {\partial u}{\partial \nu}\frac{\partial v}{\partial \nu}dS, \end{aligned} \label{e11} \begin{aligned} I_3 &= \frac 1{N_1}\int_{\Omega}\Big[(x, \nabla_x u)\Delta_y v + (x, \nabla_x v)\Delta_y u\Big]|x|^{2k}\,dx\,dy \\ &= \frac{N_1 + 2k}{N_1} \int_{\Omega} |x|^{2k} \nabla_y u\nabla_y v \,dx\,dy + \frac 1{N_1} \int_{\partial\Omega} (x, \nu_x)|x|^{2k}|\nu_y|^2 \frac{\partial u}{\partial \nu}\frac{\partial v}{\partial \nu}dS, \end{aligned} \label{e12} and \begin{aligned} I_4&=\frac{\beta}{N_2}\int_{\Omega}\Big[(y, \nabla_y u)\Delta_yv + (y, \nabla_y v)\Delta_yu\Big]|x|^{2k}\,dx\,dy\\ &= \beta\frac{N_2-2}{N_2}\int_{\Omega} |x|^{2k}\nabla_y u\nabla_y v \,dx\,dy + \frac{\beta}{N_2}\int_{\partial \Omega} |\nu_y|^2|x|^{2k}(y, \nu_y) \frac{\partial u}{\partial\nu}\frac{\partial v}{\partial\nu}dS. \end{aligned} \label{e13} Combining \eqref{e10}, \eqref{e11}, \eqref{e12} and \eqref{e13}, we obtain \begin{align*} I &= \int_{\Omega} \Big[(\frac {(x, \nabla_x u)}{N_1} + \frac{\beta(y, \nabla_y u)}{N_2})G_k v + (\frac {(x, \nabla_x v)}{N_1} + \frac{\beta(y, \nabla_y v)}{N_2})G_k u \Big]\,dx\,dy \\ &=(\beta + \frac{N_1 - 2}{N_1})\int_{\Omega} \nabla_x u \nabla_x v \,dx\,dy + (\frac {N_1 + 2k}{N_1} + \beta\frac{N_2 - 2}{N_2})\int_{\Omega} |x|^{2k} \nabla_y u\nabla_y v \,dx\,dy \\ &\quad + \frac 1{N_1} \int_{\partial \Omega} (x, \nu_x)(|\nu_x|^2 + |x|^{2k}|\nu_y|^2)\frac{\partial u}{\partial\nu}\frac{\partial v}{\partial\nu}dS\\ &\quad + \frac{\beta}{N_2}\int_{\partial\Omega}(y, \nu_y) (|\nu_x|^2+|x|^{2k} |\nu_y|^2)\frac{\partial u}{\partial \nu}\frac{\partial v}{\partial\nu}dS. \end{align*} Choosing $\beta=\frac{N_2}{N_1}(k+1)$ and taking \eqref{e5} into account, we have \begin{align*} &\frac{N(k)}{N_1}\int_{\Omega} H(x,y,u,v) \,dx\,dy\\ &= \frac{N(k) - 2}{N_1}\int_{\Omega} \Big(\nabla_x u\nabla_x v + |x|^{2k}\nabla_yu \nabla_y v\Big)\,dx\,dy \\ &\quad - \frac 1{N_1}\int_{\Omega}\Big[(x, \nabla_x H) + (k+1)(y, \nabla_y H)\Big] \,dx\,dy \\ &\quad + \frac 1{N_1}\int_{\partial\Omega}[(x, \nu_x) + (k+1)(y, \nu_y)](|\nu_x|^2 + |x|^{2k}|\nu_y|^2)\frac{\partial u}{\partial\nu}\frac{\partial v}{\partial\nu}dS. \end{align*} Noting that \begin{align*} \int_{\Omega} \Big(\nabla_x u\nabla_x v + |x|^{2k}\nabla_yu \nabla_y v\Big)\,dx\,dy &= \int_{\Omega} vf(x,y,u,v) \,dx\,dy \\ &= \int_{\Omega} ug(x,y,u,v) \,dx\,dy, \end{align*} the conclusion of Lemma \ref{lm1} follows for all $t\in \mathbb{R}$. \end{proof} The next theorem follows from Lemma \ref{lm1}. \begin{theorem} \label{thm2} Let $\Omega$ be a k-starshaped domain with respect to $\{0,0\}$ and let (S1) hold. If there exists $t\in \mathbb{R}$ such that \begin{align*} &\frac{1}{N(k)-2} [N(k)H(x, y, u, v) + (x, \nabla_x H)+(k+1)(y, \nabla_y H)]\\ &< tug(x,y,u,v) + (1-t)vf(x,y,u,v), \end{align*} for all $(x,y)\in \Omega$ and $(u,v)\in \mathbb{R}^2$, then system \eqref{e1} has no positive solution in $H^2(\Omega)\times H^2(\Omega)$. \end{theorem} In the following theorem, we get a result similar to \cite[Theorem 3.1]{g2}. \begin{theorem} \label{thm3} Let $\Omega$ be a k-starshaped domain with respect to $\{0,0\}$. If the problem \begin{gather*} -\Delta_x u - |x|^{2k}\Delta_y u = |v|^{p-1}v, \quad \text{in } \Omega, \\ -\Delta_x v - |x|^{2k}\Delta_y v = |u|^{q-1}u, \quad \text{in } \Omega, \\ u = v \quad \text{on } \partial\Omega, \end{gather*} has a nontrivial solution in $H^2(\Omega)\times H^2(\Omega)$ for $p,q \ge 1$ then $$\frac 1{p+1} + \frac 1{q+1}\ge \frac{N(k)-2}{N(k)}.$$ \end{theorem} \begin{proof} Since $\Omega$ is a k-starshaped domain, the following inequality results from Lemma \ref{lm1}, $$\frac {N(k)}{N(k)-2}\int_{\Omega} \Bigl( \frac 1{p+1} |v|^{p+1} + \frac 1{q+1} |u|^{q+1}\Bigr) \,dx\,dy \ge \int_{\Omega} \Bigl[t|v|^{p+1} + (1-t) |u|^{q+1}\Bigr]\,dx\,dy,$$ for all $t\in \mathbb{R}$. Then, we have the statement of Theorem \ref{thm3} from the fact that $$\int_{\Omega} |v|^{p+1} \,dx\,dy = \int_{\Omega} |u|^{q+1} \,dx\,dy = \int_{\Omega} (\nabla_x u\nabla_x v + |x|^{2k} \nabla_y u \nabla_y v) \,dx\,dy.$$ \end{proof} \section{Existence of classical solutions} Unlike the Laplace operator, the Grushin operator $G_k$ is not positively definite (in the domain with origin) and not radially symmetric. However, it's easy to check that, the weak maximum principle can still be applied. \begin{proposition}[Weak maximum principle] \label{prop4} Suppose that $\Omega$ is bounded. If $u \in C^2(\Omega )\cap C(\bar{\Omega})$ and $G_k u\ge 0$ in $\Omega$, then the maximum of $u$ is attained at the boundary $\partial \Omega$. \end{proposition} \begin{proof} Let $u \in C^2(\Omega )\cap C(\bar{\Omega})$. First we prove for the case $G_k u>0$ in $\Omega$ and in the next step we proceed for the case $G_k u\ge 0$ in $\Omega$. Assume $G_k u>0$ in $\Omega$ and $u$ has a maximum at $(x_0, y_0)\in \Omega$. Then \begin{gather*} \frac{\partial^2 u}{\partial x_i^2}(x_0, y_0) \le 0, \quad \text{for } i=1,\dots,N_1, \\ \frac{\partial^2 u}{\partial y_j^2}(x_0, y_0) \le 0, \quad \text{for } j=1,\dots,N_2. \end{gather*} Hence $$G_k u(x_0, y_0) = \sum_{i=1}^{N_1}\frac{\partial^2 u}{\partial x_i^2}(x_0, y_0) + |x|^{2k}\sum_{j=1}^{N_2} \frac{\partial^2 u}{\partial y_j^2}(x_0, y_0) \le 0.$$ This contradiction implies $$\sup_{\Omega} u \le \sup_{\partial\Omega} u.$$ Now we prove for the case $G_k u\ge 0$ in $\Omega$. Suppose that $\Omega\subset \{|x_1|0$ and $\alpha>0$. We have $$G_k w = G_k u + \varepsilon\alpha^2 e^{\alpha x_1} > 0 \quad \text{in } \Omega.$$ >From the arguments in first step and the construction of $w$, we deduce $$\sup_{\Omega} u \le \sup_{\Omega} w \le \sup_{\partial\Omega} w \le \sup_{\partial\Omega} u + \varepsilon e^{\alpha d}.$$ The result follows for $\varepsilon \to 0$. \end{proof} As a consequence, we have the following result. \begin{corollary} \label{coro5} If $u\in C^2(\Omega)\cap C(\bar\Omega)$ satisfies $$\begin{gathered} -G_k u \ge 0 \quad \text{in } \Omega, \\ u \ge 0 \quad \text{on } \partial\Omega, \end{gathered}\label{e14}$$ then $u\ge 0$ for $x\in \Omega$. \end{corollary} Now if $u_1, u_2 \in C^2(\Omega)\cap C(\bar\Omega)$ are solutions for $$\begin{gathered} -G_k u = f \quad \text{in } \Omega,\\ u = 0 \quad \text{on } \partial\Omega, \end{gathered} \label{e15}$$ where $f\in C(\bar\Omega)$, then $u_1-u_2$ and $u_2 - u_1$ satisfy \eqref{e14}. In other words, $u_1=u_2$ in $\Omega$. The conclusion is that \eqref{e15} has at most one solution in $C^2(\Omega)\cap C(\bar\Omega)$. \begin{proposition} \label{prop6} If $u\in C^2(\Omega)\cap C(\bar\Omega)$ is a solution of \eqref{e15} for $f\in C(\bar\Omega)$, then there exists a positive constant $C$ such that for all $(x,y)\in \Omega$, $$|u(x,y)|\le C\sup_{\Omega}|f|\,.$$ \end{proposition} \begin{proof} Denote $$\ell = \sup_{\Omega} |x|, \quad M = \sup_{\Omega} |f|, \quad v(x,y) = \frac{\ell^2-|x|^2}{2N_1}.$$ Then $G_k v = -1$ and $v\ge 0$ in $\Omega$. Put $v_1(x,y)=u(x,y)-Mv(x,y)$ and $v_2(x,y)=u(x,y)+Mv(x,y)$. By the fact that $G_k v_1=f + M\ge 0$ in $\Omega$ and using weak maximum principle, we have $$v_1(x,y) \le \sup_{\partial\Omega} v_1 \quad \text{in } \Omega,$$ or $$u(x,y) - Mv(x,y)\le \sup_{\partial\Omega}(u-Mv)\le \sup_{\partial\Omega}u=0, \quad \text{for all } (x,y)\in\Omega.$$ It follows $$u(x,y)\le Mv(x,y) \le \frac {\ell^2}{2N_1}\sup_{\Omega}|f| \quad \text{in } \Omega.$$ Arguing similarly for $v_2$ with $G_k v_2\le 0$ in $\Omega$, we get $$u(x,y)\ge \min_{\partial\Omega} u - Mv(x,y)\ge - \frac{\ell^2}{2N_1} \sup_{\Omega}|f| \quad\text{in } \Omega.$$ The proof is complete with $C=\ell^2/2N_1$. \end{proof} Before constructing the inverse operator of $G_k$, we state some conditions on the linear problem \eqref{e15}. \begin{itemize} \item[(S2)] Any solution $u$ of \eqref{e15} with $f\in C(\bar\Omega)$ belongs to $C^2(\Omega)\cap C(\bar\Omega)$. \end{itemize} \noindent {\bf Remark.} In the case $k=0$ and $\partial \Omega$ is $C^2$, hypothesis (S2) is satisfied obviously. The explicit conditions ensured the truth of (S2) for the case $k>0$ are still open. It's expected to return this problem in the other works. Given hypothesis (S2) and the uniqueness of problem \eqref{e15}, we can define the inverse operator \begin{align*} G_k^{-1}: C(\bar\Omega) & \to C(\bar\Omega), \\ f & \mapsto u. \end{align*} Furthermore, Proposition \ref{prop6} ensures that $G_k^{-1}$ is compact. \begin{itemize} \item[(S3)] The functions $f(x,y,s,t)$, $g(x,y,s,t)$ in $C(\bar\Omega\times \mathbb{R} \times \mathbb{R})$ are nondecreasing in $s$ and $t$ for all $(x,y)\in\Omega$, i.e. the maps $s\mapsto f(x,y,s,t), s\mapsto g(x,y,s,t)$ for fixed $t\in \mathbb{R}$ and $t\mapsto f(x,y,s,t), t\mapsto g(x,y,s,t)$ for fixed $s\in \mathbb{R}$ are nondecreasing for all $(x,y)\in \Omega$. \end{itemize} \noindent {\bf Definition.} %2 A pair $(\bar u, \bar v)\in C^2(\Omega)\cap C(\bar\Omega)$ is said to be a supersolution to \eqref{e1} if $$\begin{gathered} - G_k \bar u \ge f(x,y,\bar u,\bar v) \quad\text{in }\Omega,\\ - G_k \bar v \ge g(x,y,\bar u,\bar v) \quad\text{in }\Omega,\\ \bar u \ge 0,\quad \bar v \ge 0 \quad\text{on }\partial\Omega. \end{gathered} \label{e17}$$ Similarly, $(\underline u, \underline v)\in C^2(\Omega)\cap C(\bar\Omega)$ is called subsolution to \eqref{e1} if $$\begin{gathered} - G_k \underline u \le f(x,y,\underline u,\underline v) \quad\text{in }\Omega,\\ - G_k \underline v \le g(x,y,\underline u,\underline v) \quad\text{in }\Omega,\\ \underline u \le 0, \quad \underline v \le 0 \quad\text{on }\partial\Omega. \end{gathered} \label{e18}$$ The following theorem is the main result for current section. \begin{theorem} \label{thm7} Assume that hypotheses (S2) and (S3) hold. If \eqref{e1} has a subsolution ($\underline u, \underline v$) and a supersolution ($\bar u, \bar v$) such that $(\underline u, \underline v) \le (\bar u, \bar v)$, then there exists at least one solution ($u, v$) of \eqref{e1} satisfying $(\underline u, \underline v) \le (u, v) \le (\bar u, \bar v)$. \end{theorem} \begin{proof} We use similar arguments as in \cite{d1}. Let $f^{*}(x, y, u, v)$ and $g^{*}(x,y,u,v)$ be the functions which are defined from $f(x, y, u,v)$, $g(x,y,u,v)$ as follows \begin{itemize} \item[(i)] If $u<\underline u$ then replace $u$ by $\underline u$, if $u>\bar u$ then replace $u$ by $\bar u$. \item[(ii)] If $v<\underline v$ then replace $v$ by $\underline v$, if $v>\bar v$ then replace $v$ by $\bar v$. \end{itemize} Then $f^{*}$, $g^{*}$ are continuous functions. Consider the problem $$\begin{gathered} - G_k u = f^{*}(x,y,u,v) \quad\text{in }\Omega, \\ - G_k v = g^{*}(x,y,u,v) \quad\text{in }\Omega, \\ u = v = 0 \quad\text{on }\partial\Omega. \end{gathered}\label{e19}$$ The definition of $f^{*}(x, y, u, v)$ and $g^{*}(x,y,u,v)$ implies that, if $(u,v)$ is a solution of \eqref{e19} satisfying $$(\underline u, \underline v)\le (u,v) \le (\bar u, \bar v), \label{e20}$$ then $(u, v)$ is a solution for \eqref{e1}. Let us define the operator $$T(u,v) = \begin{pmatrix} -G_k^{-1} & 0 \\ 0 & -G_k^{-1} \end{pmatrix} \begin{pmatrix} f^{*}(x,y,u,v) \\ g^{*}(x,y,u,v) \end{pmatrix}$$ on the convex bounded subset of $[C(\overline\Omega)]^{2}$ $$D = \{(u,v): (\underline u, \underline v)\le (u,v)\le (\bar u, \bar v)\}.$$ Since $f^{*}$ and $g^{*}$ are continuous and bounded, $G_k^{-1}$ is compact, hence $T$ is compact. Note now that system \eqref{e19} is equivalent to $(u,v)=T(u,v)$. To apply Shauder fixed point theorem \cite[page 60]{g2}, we have only to prove that $T(D)\subset D$. Indeed, let $(u,v)$ be a solution to \eqref{e19}, we show that $(u,v)$ satisfies \eqref{e20}. It suffices to prove that, the following sets \begin{gather*} \Omega_1 = \{(x, y)\in \Omega: u(x,y) > \bar u(x,y)\}, \\ \Omega_2 = \{(x, y)\in \Omega: v(x,y) > \bar v(x,y)\}, \\ \Omega_3 = \{(x, y)\in \Omega: u(x,y) < \underline u(x,y)\}, \\ \Omega_4 = \{(x, y)\in \Omega: v(x,y) < \underline v(x,y)\} \end{gather*} are empty. We proceed with $\Omega_1$, the proofs for the other set are performed analogously. Assume to the contrary that $\Omega_1$ is not empty. The continuity of $u$ and $\bar u$ ensures that $\Omega_1$ is open set in $\Omega$, then $\Omega_1 \cap \partial\Omega_1 = \emptyset$. The first inequality in \eqref{e17} and the first equation in \eqref{e19} gives $$- G_k (u-\bar u)\le f^{*}(x, y, u,v) - f(x,y,\bar u,\bar v). \label{e21}$$ For pointing out the sign of the left hand side in \eqref{e21}, we define the following subsets of $\Omega_1$ \begin{gather*} \Omega_1^{+} = \{(x,y)\in \Omega_1: v(x,y) > \bar v(x,y)\}, \\ \Omega_1^{-} = \{(x,y)\in \Omega_1: v(x,y) < \underline v(x,y)\}, \\ \Omega_1^{0} = \{(x,y)\in \Omega_1: \underline v(x,y) \le v(x,y) \le \bar v(x,y)\}. \end{gather*} It's easy to see that $\Omega_1 = \Omega_1^{+}\cup \Omega_1^{-}\cup \Omega_1^{0}$ and each subset is separated from others. By the definition of $f^{*}$ and monotoncity of $f$, we have \begin{gather*} f^{*}(x, y, u,v) - f(x,y,\bar u,\bar v) = f(x,y,\bar u, \bar v) - f(x,y,\bar u, \bar v)=0 \quad\text{in }\Omega_1^{+}, \\ f^{*}(x, y, u,v) - f(x,y,\bar u,\bar v) = f(x,y,\bar u, \underline v) - f(x,y,\bar u, \bar v) \le 0 \quad\text{in }\Omega_1^{-}, \\ f^{*}(x, y, u,v) - f(x,y,\bar u,\bar v) = f(x,y,\bar u, v) - f(x,y,\bar u, \bar v)\le 0 \quad\text{in }\Omega_1^{0}. \end{gather*} Taking \eqref{e21} into account, we conclude that $G_k (u-\bar u)\ge 0$. Proposition \ref{prop4} concludes that, the maximum of $u-\bar u$ is attained at $(x^1, y^1)\in \partial\Omega_1$. Therefore, $(u-\bar u)(x^1, y^1)>0$. The contradiction occurs since $(x_1, y_1)\in \Omega_1 \cap \partial\Omega_1$. \end{proof} \subsection*{The critical cases} We construct a class of subsolutions and supersolutions of \eqref{e1} in the case when nonlinearities have critical growth $$\begin{gathered} \Delta_x u + |x|^{2k}[\Delta_y u + |v|^{\frac{N(k)+2}{N(k) - 2}} + h_1(x,y)] = 0 \quad\text{in }\Omega, \\ \Delta_x v + |x|^{2k}[\Delta_y v + |u|^{\frac{N(k)+2}{N(k) - 2}} + h_2(x,y)] = 0 \quad\text{in }\Omega, \\ u = v = 0 \quad\text{on }\partial\Omega, \end{gathered}\label{e22}$$ where $h_1, h_2 \in C(\bar\Omega)$. Denote $$\rho = (|x|^{2k+2} + (k+1)^2|y|^2)^{\frac 1{2k+2}}.$$ We find a class of supersolutions $(u^{*}, v^{*})$ to \eqref{e22} of the form $$u^{*} = v^{*} = C\lambda^{\frac{N(k)-2}2} (1 + \lambda^2\rho^2)^{-\frac{N(k)-2}2}, \label{e23}$$ where $C>0$, $\lambda > 0$. Note that, if $u=u(\rho)$ then $$\Delta_x u + |x|^{2k}\Delta_y u = \Big[\frac{\partial^2 u}{\partial \rho^2} + \frac{N(k)-1}{\rho}\frac{\partial u}{\partial\rho}\Big]\frac{|x|^{2k}}{\rho^{2k}}.$$ It suffices to find $C$ such that \begin{gather*} \frac{\partial^2 u^*}{\partial \rho^2} + \frac{N(k)-1}{\rho} \frac{\partial u^*}{\partial\rho} + \rho^{2k}(|v^*|^{\frac{N(k)+2}{N(k)-2}} + h_1(x,y)) \le 0 \quad\text{in }\Omega, \\ \frac{\partial^2 v^*}{\partial \rho^2} + \frac{N(k)-1}{\rho} \frac{\partial v^*}{\partial\rho} + \rho^{2k}(|u^*|^{\frac{N(k)+2}{N(k)-2}} + h_2(x,y)) \le 0 \quad\text{in }\Omega, \\ u \ge 0, \quad v \ge 0 \quad\text{on }\partial\Omega. \end{gather*} Put $$Au=-\Big[\frac{\partial^2 u}{\partial \rho^2} + \frac{N(k)-1}{\rho} \frac{\partial u}{\partial\rho} + \rho^{2k}|u|^{\frac{N(k)+2}{N(k)-2}}\Big].$$ From \eqref{e23} and with some computations (Using Maple software), we have $$A u^* = -C\lambda^{\frac{N(k)+2}2}(1+\lambda^2 \rho^2)^{-\frac{N(k)+2}2} [C^{\frac 4{N(k)-2}}\rho^{2k} - N(k)(N(k)-2)].$$ Choosing $$C< [N(k)(N(k)-2)]^{\frac {N(k)-2}4}(\max_{\bar\Omega}\rho^{2k})^{-\frac{N(k)-2}4},$$ we have $A u^* > 0$ in $\Omega$. Obviously, $(u^*, v^*)$ is a supersolution and $(0, 0)$ is a subsolution to \eqref{e22} if $$0 < h_1(x,y), h_2(x,y) \le (\max_{\bar\Omega} \rho)^{-2k} Au^* \quad\text{in } \Omega$$ (in this case, $(0,0)$ can not be solution of \eqref{e22}). As a consequence of Theorem \ref{thm7}, there exists a positive solution $(u,v)$ for \eqref{e22} such that $0< u$, $v \le u^{*}$ for all $(x,y)\in \Omega$. \subsection*{The supercritical cases} Consider the problem $$\begin{gathered} \Delta_x u + |x|^{2k}[\Delta_y u + \alpha |v|^{p} + h_1(x,y)] = 0 \quad\text{in }\Omega, \\ \Delta_x v + |x|^{2k}[\Delta_y v + \beta |u|^{q} + h_2(x,y)] = 0 \quad\text{in }\Omega, \\ u = v = 0 \quad\text{on }\partial\Omega, \end{gathered}\label{e24}$$ where $p, q> \frac{N(k)+2}{N(k)-2}$; $\alpha, \beta > 0$; $h_1, h_2 \in C(\bar\Omega)$. One finds a class of supersolutions ($u^{\#}, v^{\#}$) to \eqref{e24} of the form \begin{gather*} u^{\#} = C_1\lambda^{a}(1+\lambda^2 \rho^2)^{-a}, \\ v^{\#} = C_2\lambda^{b}(1+\lambda^2 \rho^2)^{-b}, \end{gather*} where $C_1>0$, $C_2>0$; $\lambda > 0$; $a, b>0$. For $\gamma>0$ and $s>1$, let $$B_{s, \gamma}(u, v) = \frac{\partial^2 u}{\partial \rho^2} + \frac{N(k)-1}{\rho} \frac{\partial u}{\partial\rho} + \gamma \rho^{2k}|v|^s\,.$$ After some calculations (using Maple software) we get \begin{align*} B_{p, \alpha}(u^{\#}, v^{\#}) &= -2C_1 a\lambda^{2+a}(1+\lambda^2\rho^2)^{-2-a}[\lambda^2\rho^2(N(k)-2a-2)+N(k)]\\ &\quad + \alpha\rho^{2k} C_2^p\lambda^{bp}(1+\lambda^2\rho^2)^{-bp}. \end{align*} Similarly, \begin{align*} B_{q, \beta}(v^{\#}, u^{\#}) &= -2C_2 b\lambda^{2+b}(1+\lambda^2\rho^2)^{-2-b}[\lambda^2\rho^2(N(k)-2b-2)+N(k)]\\ &\quad + \beta\rho^{2k} C_1^q\lambda^{aq}(1+\lambda^2\rho^2)^{-aq}. \end{align*} We choose $a, b$ such that $qa = 2+b$ and $pb = 2+a$. Clearly $$a = \frac{2p+2}{pq-1}, \quad b = \frac{2q+2}{pq-1}.$$ Therefore, \begin{align*} &B_{p, \alpha}(u^{\#}, v^{\#})\\ & = \lambda^{2+a}(1+\lambda^2\rho^2)^{-2-a}\{\alpha \rho^{2k} C_2^{p} - 2aC_1 [\lambda^2\rho^2 (N(k)-2a-2) + N(k)] \}, \\ &B_{q, \beta}(v^{\#}, u^{\#}) \\ & =\lambda^{2+b}(1+\lambda^2\rho^2)^{-2-b}\{\beta \rho^{2k} C_1^{q} - 2bC_2 [\lambda^2\rho^2 (N(k)-2b-2) + N(k)] \}. \end{align*} Under the above mentioned conditions for $p$ and $q$, we have $$N(k)-2a-2 > 0, \ N(k)-2b-2 >0.$$ For the positivity of $-B_{p, \alpha}(u^{\#}, v^{\#})$ and $-B_{q, \beta}(v^{\#}, u^{\#})$, the following conditions are sufficient: \begin{gather*} \alpha \max_{\bar\Omega}\rho^{2k} C_2^p < 2 a N(k)C_1, \\ \beta \max_{\bar\Omega}\rho^{2k} C_1^q < 2 b N(k)C_2. \end{gather*} That implies \begin{gather*} C_1 < \Big[\frac{2bN(k)}{\beta \max_{\bar\Omega}\rho^{2k}}\Big]^{\frac p{pq-1}} \Big[\frac{2aN(k)}{\alpha \max_{\bar\Omega}\rho^{2k}}\Big]^{\frac 1{pq-1}}, \\ C_2 < \Big[\frac{2bN(k)}{\beta \max_{\bar\Omega}\rho^{2k}}\Big]^{\frac 1{pq-1}} \Big[\frac{2aN(k)}{\alpha \max_{\bar\Omega}\rho^{2k}}\Big]^{\frac q{pq-1}}. \end{gather*} Now, $(u^{\#}, v^{\#})$ becomes a supersolution and $(0, 0)$ becomes a subsolution to \eqref{e24} if \begin{gather*} 0< h_1(x,y) \le - B_{p,\alpha}(u^{\#}, v^{\#})(\max_{\bar\Omega}\rho)^{-2k}, \\ 0< h_2(x,y) \le - B_{q,\beta}(v^{\#}, u^{\#})(\max_{\bar\Omega}\rho)^{-2k}. \end{gather*} Applying Theorem \ref{thm7}, we see that there exists at least one solution $(u, v)$ for \eqref{e24} satisfying \begin{gather*} 0 < u \le u^{\#}, \\ 0 < v \le v^{\#}. \end{gather*} \section{Maximal and minimal weak solutions} \noindent {\bf Definition.} %3 By $S_1^p(\Omega)$, $1 0$ and $s\ge 0$, provided small positive number $\tau$. So, one can take the norm for $S_{1,0}^2(\Omega)$ as follows $$\|(u,v)\|_{S_{1,0}^2(\Omega)}^2 = \int_{\Omega}\Big[|\nabla_x u|^2 + |\nabla_x v|^2 + |x|^{2k}(|\nabla_y u|^2 + |\nabla_y v|^2)\Big] \,dx\,dy.$$ The definition of the weak solution for system \eqref{e1} is stated as follows: \noindent {\bf Definition.} A pair of functions $(u,v)\in S^2_{1,0}(\Omega)$ is called weak solution for system \eqref{e1} if \begin{gather*} \int_{\Omega} [\nabla_x u \nabla_x\varphi + |x|^{2k}\nabla_y u\nabla_y\varphi] \,dx\,dy = \int_{\Omega} f(x,y,u,v)\varphi \,dx\,dy, \\ \int_{\Omega} [\nabla_x v \nabla_x\psi + |x|^{2k}\nabla_y v\nabla_y\psi]\,dx\,dy = \int_{\Omega} g(x,y,u,v)\psi \,dx\,dy \end{gather*} for all $\varphi, \psi \in C^1_{0}(\Omega)$. \smallskip The following definition describes the comparison on boundary $\partial \Omega$ of two functions in $S^2_1(\Omega)$. \noindent {\bf Definition} %5 Let $(u, v) \in S^2_1(\Omega)$. A function $v$ is said to be less than or equal a function $u$ on $\partial\Omega$ if $(max\{0, v-u\}, 0) \in S^2_{1, 0}(\Omega)$. \smallskip The proof of the following assertion is standard and therefore omitted. \begin{proposition} \label{prop8} Consider the pair of functions $(u_1, u_2)\in S_{1}^2(\Omega)\cap (L^{\infty}(\Omega))^2$ such that, for all $\varphi$ satisfying $\varphi \in C^1_0 (\Omega), \varphi \ge 0$ in $\Omega$, \begin{gather*} \int_{\Omega} [\nabla_x u_1 \nabla_x \varphi + |x|^{2k} \nabla_y u_1\nabla_y \varphi]\,dx\,dy \le \int_{\Omega} [\nabla_x u_2 \nabla_x \varphi + |x|^{2k} \nabla_y u_2 \nabla_y \varphi]\,dx\,dy, \\ u_1 \le u_2 \quad\text{on }\partial\Omega. \end{gather*} Then $u_1\le u_2$ a.e in $\Omega$. \end{proposition} \begin{proposition} \label{prop9} For every $h_1, h_2 \in L^2(\Omega)$, the problem $$\begin{gathered} -G_k u = h_1(x,y), \quad (x,y) \in \Omega, \\ -G_k v = h_2(x,y), \quad (x,y) \in \Omega, \\ u = v = 0, \quad (x,y) \in \partial\Omega, \end{gathered}\label{e25}$$ admits a unique weak solution $(u, v) \in S^2_{1,0}(\Omega)$. Moreover, the associated operator $W: (h_1, h_2) \mapsto (u, v)$ is continuous and nondecreasing. \end{proposition} \begin{proof} The problem \eqref{e25} is written in the form of system in order to construct operator $W$. The proof of existence and uniqueness can be proceeded using the Riesz representation theorem in Hilbert space $S^2_{1,0}(\Omega)$ (as in \cite[page 150]{c3}). The continuous property of $W$ follows from the estimate $$\|(u,v)\|_{S^2_{1,0}(\Omega)}\le C\|(h_1, h_2)\|_{(L^2(\Omega))^2} = C\Big[\int_{\Omega}(|h_1(x,y)|^2 + |h_2(x,y)|^2)\,dx\,dy\Big]^{1/2},$$ where $C$ is a positive constant. Assume that $(h_1, h_2), (l_1, l_2) \in (L^2(\Omega))^2$ and $(u_h, v_h)$, $(u_l, v_l)\in S^2_{1,0}(\Omega)$ are solutions of \eqref{e25} respectively, then we have \begin{gather*} \int_{\Omega} \Big[(l_1-h_1)\varphi \Big]\,dx\,dy = \int_{\Omega}\Big[\nabla_x (u_l-u_h) \nabla_x \varphi + |x|^{2k} \nabla_y (u_l-u_h)\nabla_y \varphi\Big]\,dx\,dy, \\ \int_{\Omega} \Big[(l_2-h_2)\psi \Big]\,dx\,dy = \int_{\Omega}\Big[\nabla_x (v_l-v_h) \nabla_x \psi + |x|^{2k} \nabla_y (v_l-v_h) \nabla_y \psi\Big]\,dx\,dy, \end{gather*} for all $\varphi, \psi \in C^1_0 (\Omega)$ satisfying $\varphi, \psi\ge 0$ in $\Omega$. Applying Proposition \ref{prop8}, we obtain the nondecreasing property of $W$. \end{proof} Before stating the results for this section, we replace hypothesis (S3) by \begin{itemize} \item[(S3')] $f(x,y,s,t), g(x,y,s,t)$ are Caratheodory functions: $f(x,y,.,.), g(x,y,.,.)$ are continuous for a.e. $(x,y)\in \Omega$, $f(.,.,s,t), g(.,.,s,t)$ are measurable for all $(s,t) \in\mathbb{R}^2$ and $f(.,.,s,t), g(.,.,s,t)$ are bounded if $(s,t)$ belong to bounded sets. In addition, $f(x,y,.,t), g(x,y,.,t)$ for fixed $t\in\mathbb{R}$ and $f(x,y,s,.), g(x,y,s,.)$ for fixed $s\in \mathbb{R}$ are nondecreasing for a.e. $(x,y) \in \Omega$. \end{itemize} Now let us define the subsolutions and supersolutions for \eqref{e1} in the weak sense. The comparison on $\partial \Omega$ is realized by the definition above. \smallskip \noindent {\bf Definition.} %6 Let $(\underline u, \underline v), (\overline u, \overline v) \in S^2_{1}(\Omega) \cap (L^{\infty}(\Omega))^2$. These pairs of functions are said to be a system of subsolution and supersolution in the weak sense for \eqref{e1} if they satisfy: \begin{itemize} \item[(a)] $\underline u(x,y) \le \overline u(x,y)$, $\underline u(x,y) \le \overline u(x,y)$ a.e. in $\Omega$, $\underline u(x,y) \le 0\le \overline u(x,y)$, $\underline u(x,y) \le 0\le \overline u(x,y)$ on $\partial\Omega$, \item[(b)] For all $\varphi\in C^1_0(\Omega): \varphi \ge 0 \quad\text{in }\Omega$, \begin{gather*} \int_{\Omega} [\nabla_x \overline u\nabla_x\varphi + |x|^{2k}\nabla_y \overline u\nabla\varphi]\,dx\,dy \ge \int_{\Omega} f(x,y,\overline u, \overline v)\varphi \,dx\,dy,\\ \int_{\Omega} [\nabla_x \underline u\nabla_x\varphi + |x|^{2k}\nabla_y \underline u\nabla\varphi]\,dx\,dy \le \int_{\Omega} f(x,y,\underline u, \underline v)\varphi \,dx\,dy, \end{gather*} and \begin{gather*} \int_{\Omega} [\nabla_x \overline v\nabla_x\varphi + |x|^{2k}\nabla_y \overline v\nabla\varphi]\,dx\,dy \ge \int_{\Omega} g(x,y, \overline u,\overline v)\varphi \,dx\,dy,\\ \int_{\Omega} [\nabla_x \underline v\nabla_x\varphi + |x|^{2k}\nabla_y \underline v\nabla\varphi]\,dx\,dy \le \int_{\Omega} g(x,y, \underline u, \underline v)\varphi \,dx\,dy. \end{gather*} \end{itemize} \smallskip The following theorem is the main result of this section. \begin{theorem} \label{thm10} Assume (S3') and let $(\underline u, \underline v), (\overline u, \overline v)$ be a system of subsolution and supersolution of \eqref{e1}. Then, there exists a minimal (and, respectively, a maximal) weak solution $(u_*, v_*)$ (respectively $(u^*, v^*)$) for problem \eqref{e1} in the set \begin{align*} [(\underline u, \underline v), (\overline u, \overline v)] = \big\{&(u, v)\in (L^{\infty}(\Omega))^2: \underline u(x,y) \le u(x, y) \le \overline u(x,y), \\ &\underline v(x,y) \le v(x, y)\le \overline v(x,y), \text{ a.e. in } \Omega \big\}. \end{align*} Precisely, every weak solution $(u,v)\in [(\underline u, \underline v), (\overline u, \overline v)]$ of \eqref{e1} satisfies $u_*(x,y)\le u(x,y) \le u^*(x,y), v_*(x,y)\le v(x,y)\le v^*(x,y)$ for a.e. $(x, y)\in \Omega$. \end{theorem} \begin{proof} Note that the set $[(\underline u, \underline v), (\overline u, \overline v)]$ is a subset of $(L^{\infty}(\Omega))^2$ with the topology of a.e. convergence. We define the operator $Z: [(\underline u, \underline v), (\overline u, \overline v)] \to (L^2(\Omega))^2$ by $$Z(\tilde u, \tilde v) = \Big( f(.,., \tilde u(.,.), \tilde v(.,.)), \ g(.,., \tilde u(.,.), \tilde v(.,.))\Big).\label{e26}$$ By hypothesis (S3'), we get that $Z$ is nondecreasing and bounded. Moreover, if $(\tilde u_n, \tilde v_n), (\tilde u, \tilde v)$ is in $[(\underline u, \underline v), (\overline u, \overline v)]$ then \begin{align*} \|Z(\tilde u_n, \tilde v_n) - Z (\tilde u, \tilde v) \|_{(L^2(\Omega))^2}^2 &= \int_{\Omega} |f(x,y,\tilde u_n, \tilde v_n) - f(x,y,\tilde u, \tilde v)|^2 \,dx\,dy\\ &\quad + \int_{\Omega} |g(x,y,\tilde u_n, \tilde v_n) - g(x,y,\tilde u, \tilde v)|^2 \,dx\,dy. \end{align*} Let $(\tilde u_n, \tilde v_n) \to (\tilde u, \tilde v)$ a.e in $\Omega$. Applying the Lebesgue dominated theorem, we obtain that $\|Z(\tilde u_n, \tilde v_n) - Z (\tilde u, \tilde v) \|_{(L^2(\Omega))^2} \to 0$ and $Z$ is continuous. The constructions of the operator $W$ in Proposition \ref{prop9} and operator $Z$ in \eqref{e26} allow us to define the operator $T: [(\underline u, \underline v), (\overline u, \overline v)]\to S^2_{1,0}(\Omega)$ by $T=W\circ Z$. It's easy to see that, for a pair $(\tilde u, \tilde v)$ in $[(\underline u, \underline v), (\overline u, \overline v)], T(\tilde u, \tilde v)$ is the unique weak solution of the boundary-value problem $$\begin{gathered} -G_k u = f(x,y,\tilde u, \tilde v), \quad (x,y) \in \Omega, \\ -G_k v = g(x,y,\tilde u, \tilde v), \quad (x,y) \in \Omega, \\ u = v = 0, \quad (x,y) \in \partial\Omega. \end{gathered}\label{e27}$$ Putting $(u_1, v_1)=T(\underline u, \underline v)$, $(u^1, v^1)=T(\overline u, \overline v)$, we deduce that, for all $\varphi$ satisfying $\varphi \in C^1_0(\Omega)$, $\varphi\ge 0$ in $\Omega$, \begin{align*} \int_{\Omega} [\nabla_x u_1\nabla_x\varphi + |x|^{2k}\nabla_y u_1\nabla_y \varphi]\,dx\,dy &= \int_{\Omega} f(x,y,\underline u, \underline v)\varphi \,dx\,dy\\ &\ge \int_{\Omega} [\nabla_x \underline u \nabla_x\varphi + |x|^{2k}\nabla_y \underline u\nabla_y \varphi]\,dx\,dy, \end{align*} \begin{align*} \int_{\Omega} [\nabla_x v_1\nabla_x\varphi + |x|^{2k}\nabla_y v_1\nabla_y \varphi]\,dx\,dy & = \int_{\Omega} g(x,y,\underline u, \underline v)\varphi \,dx\,dy\\ &\ge \int_{\Omega} [\nabla_x \underline v \nabla_x\varphi + |x|^{2k}\nabla_y \underline v\nabla_y \varphi]\,dx\,dy. \end{align*} By applying Proposition \ref{prop8}, we obtain $\underline u\le u_1, \underline v\le v_1$, or briefly, $(\underline u, \underline v)\le T(\underline u,\underline v)$. By the same arguments, we get that $T(\overline u, \overline v)\le (\overline u, \overline v)$. Taking into account that $T$ is nondecreasing, we have $$(\underline u, \underline v)\le T(\underline u,\underline v)\le T(u,v)\le T(\overline u, \overline v)\le (\overline u, \overline v),$$ for all $(u,v)\in [(\underline u, \underline v), (\overline u, \overline v)]$. It's now to construct two sequences $(u_n, v_n)$ and $(u^n, v^n)$ as follows \begin{gather*} (u_0, v_0) =(\underline u, \underline v), \ (u_{n+1}, v_{n+1}) = T(u_n, v_n), \\ (u^0, v^0) =(\overline u, \overline v), \ (u^{n+1}, v^{n+1}) = T(u^n, v^n). \end{gather*} Repeating the same process, we can prove that $$(u_0, v_0)\le (u_1, v_1)\le\dots\le(u_n, v_n)\le (u,v) \le (u^n, v^n)\le \dots\le(u^1, v^1)\le (u^0, v^0)$$ a.e. in $\Omega$, for every weak solution $(u,v)\in [(\underline u, \underline v), (\overline u, \overline v)]$. Then $(u_n, v_n) \to (u_*, v_*)$, $(u^n, v^n)\to (u^*, v^*)$, a.e. in $\Omega$. Obviously, $(u_*, v_*)$ and $(u^*,v^*)\in (L^{\infty}(\Omega))^2$ and $(u_*, v_*), (u^*,v^*)\in [(\underline u, \underline v), (\overline u, \overline v)]$. Then, by the fact of $T(\tilde u, \tilde v)$ commented in \eqref{e27}, we have that $T(u_*,v_*)$ (respectively $T(u^*,v^*)$) is the unique weak solution of \eqref{e27} when $(\tilde u, \tilde v)$ is replaced by $(u_*,v_*)$ (respectively by $(u^*,v^*)$). Considering \eqref{e27} as the linear problem in Proposition \ref{prop9}, we have the conclusion that $T(u_*,v_*)$ and $T(u^*,v^*)\in S^2_{1,0}(\Omega)$. Since the continuity of $Z$ and $W$ ensures the continuity of $T$, we deduce that $(u_{n+1}, v_{n+1}) = T(u_n, v_n) \to T(u_*, v_*) = (u_*, v_*)$, $(u^{n+1}, v^{n+1}) = T(u^n, v^n)\to T(u^*, v^*) = (u^*, v^*)$ in $S^2_{1,0}(\Omega)$. This completes the proof. \end{proof} \subsection*{Acknowledgement} The authors would like to express their gratitude to the anonymous referee for his/her valuable suggestions and comments. \begin{thebibliography}{00} \bibitem{a1} C. O. Alves, D. C.de Morais Filho \& M. A. S. Souto, \emph{On system of Elliptic equations involving subcritical or critical Sobolev exponents}, Nonlinear Analysis, Vol. ?? (1998), 1-13. \bibitem{b1} Ahmed Bensedik and Mohamed Bouchekif, \emph{On certain nonlinear elliptic systems with indefinite terms}, Electron. J. Diff. Eqns., Vol 2002 (2002) No. 83, 1-16. \bibitem{c1} N. M. Chuong, \emph{Degenerate parabolic pseudodifferential operators of variable order}, Dokl. Akad. Nauk SSSR, Vol. 268 (1983) No. 5, 1055-1058. \bibitem{c2} N. M. Chuong, T.D. Ke, N. V. Thanh, and N. M. Tri, \emph{Non-existence theorems for boundary value problems for some classes of semilinear degenerate elliptic operators}, Proceedings of the conference on Partial Differential Equations and their Applications, Hanoi Dec. 1999, (2001), 185-190. \bibitem{c3} N. M. Chuong, H. T. Ngoan, N. M. Tri, and L. Q. Trung, \emph{Partial Differential Equations} Education Publishing House (2000), MR 2001i:35001. \bibitem{d1} E. N. Dancer and G. Sweers, \emph{On the existence of maximal weak solution for a semilinear elliptic equation}, Differential and Integral Equations, Vol. 2 (1989) No. 4, 533-540. \bibitem{d2} Klaus Deimling, \emph{Nonlinear Functional Analysis} Springer-Verlag Berlin Heidelberg, (1985). \bibitem{g1} V. V. Grushin, \emph{On a class of Elliptic Pseudo Differential Operators Degenerate on a Submanifold}, Math. USSR Sbornik, Vol. 13 (1971), 155-183. \bibitem{g2} M. Garcia-Huidobro, R. Manasevich, E. Mitidieri and C. Yarur, \emph{Existence and Nonexistence of Positive Singular Solutions for a Class of Semilinear Elliptic Systems}, Arch. Rational Mech. Anal., Vol. 140 (1997), 253-284. \bibitem{k1} T. D. Ke, \emph{Existence of non-negative solutions for a semilinear degenerate elliptic system} Proceedings of the international conference on Abstract and Applied Analysis (edited by N. M. Chuong, L. Nirenberg, L. H. Son, W. Tutschke), Hanoi Aug. 2002 (World Scientific 2004) (to appear). \bibitem{t1} N. T. C. Thuy and N. M. Tri, \emph{Existence and Nonexistence Results for Boundary Value Problems (BVP) for Semilinear Elliptic Degenerate Operators}, Russian Journal of Mathematical Physics, Vol. 9 (2002), No. 3, 366-371. \bibitem{t2} Nguyen Minh Tri, \emph{On Grushin's Equation}, Matemachicheskie Zametki, Vol. 63 (1988), 95-105. \end{thebibliography} \end{document}