\documentclass[reqno]{amsart}
\usepackage{hyperref}

\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2016 (2016), No. 177, pp. 1--12.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu}
\thanks{\copyright 2016 Texas State University.}
\vspace{8mm}}

\begin{document}
\title[\hfilneg EJDE-2016/177\hfil Biharmonic systems with
 singular nonlinearity]
{Existence of solutions to biharmonic systems with
 singular nonlinearity}

\author[A. L. A. de Araujo, L. F. O. Faria \hfil EJDE-2016/177\hfilneg]
{Anderson L. A. de Araujo, Luiz F. O. Faria}

\address{Anderson L. A. de Araujo \newline
Departamento de Matem\'{a}tica,
Universidade Federal de Vi\c{c}osa,
CEP 36570-900 Vi\c{c}osa,
Minas Gerais, Brazil}
\email{anderson.araujo@ufv.br}

\address{Luiz F. O. Faria \newline
Departamento de Matem\'{a}tica - ICE,
Universidade Federal de Juiz de Fora,
CEP 36036-330 Juiz de Fora,
 Minas Gerais, Brazil}
\email{luiz.faria@ufjf.edu.br}

\thanks{Submitted March 22, 2016. Published July 6, 2016.}
\subjclass[2010]{35B09, 35B45, 35J58, 35J75}
\keywords{Biharmonic operator; singular nonlinearity; Green function; 
\hfill\break\indent Schauder's fixed point theorem}

\begin{abstract}
 In this article we prove the existence of positive solutions
 of nonlinear singular biharmonic elliptic systems in smooth bounded domains,
 with coupling of the equations, under Navier boundary condition.
 Under some suitable assumptions on the nonlinearity, we prove a uniqueness
 result. The existence result is based on the Schauder's fixed point theorem.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{remark}[theorem]{Remark}
\allowdisplaybreaks


\section{Introduction}

In this article, we study the existence and uniqueness of positive 
solutions to the biharmonic elliptic system
\begin{equation}\label{P}
\begin{gathered}
 \Delta^2u = \frac{A(x)}{u^{\alpha}} + \frac{B(x)}{(u+v)^{r_1}} 
 \quad\text{in }\Omega,\\
  \Delta^2v = \frac{C(x)}{v^{\beta}} + \frac{D(x)}{(u+v)^{r_2}} 
\quad\text{in }\Omega,\\
 u, v>0 \quad\text{in }\Omega,\\
 u=0, \quad\Delta u=0 \quad\text{on } \partial\Omega,\\
 v=0, \quad \Delta v=0 \quad\text{on } \partial\Omega,
\end{gathered}
\end{equation}
where $\Delta^2$ is the biharmonic operator, $\alpha,\beta, r_1,r_2$ 
are positive constants, $\Omega\subset\mathbb{R}^N$ ($1\leq N$)  is a 
smooth bounded domain and $A,\,B,\,C,\,D\,\in C(\overline{\Omega})$. 
The condition on the boundary is known as Navier boundary condition.

System \eqref{P} appears as a natural extension of the single singular problem
\begin{equation}\label{pa}
\begin{gathered}
 \Delta^2u = \frac{A(x)}{u^{\alpha}}  \quad\text{in }\Omega,\\
  u=0, \quad \mathfrak{B}(u)=0 \quad\text{on } \partial\Omega,\\
\end{gathered}
\end{equation}
where $0<\alpha<1$, which has been considered, among other works, 
in \cite{HC} (when $\mathfrak{B}(u)=\Delta u$). The problem \eqref{pa}
was also studied under Dirichlet boundary condition (that is, when 
$\mathfrak{B}(u)=\partial_{\nu} u$), see \cite{Ghergu}. In both references, 
the existence result was obtained by means of Schauder fixed point theorem. 
The study of singular elliptic problems is greatly justified in view of some 
basic aspects of mathematical research. They arise in several physical 
situations such as fluids, biological pattern formation and so on. 
As a physical illustration we describe briefly a practical problem which 
leads to a singular problem as it is has been done in Fulks-Maybe \cite{FM}.

The single fourth order elliptic equations arises in the study of traveling 
waves in suspension bridges \cite{LazerMcKenna}. In recent years, fourth order 
nonlinear differential equations have become increasingly popular due to their 
possible applications in the fields of image and signal processing, 
nuclear physics, and engineering, see e.g. \cite{cmm,CX,yk}.
 The current knowledge of fourth order elliptic equations has
considerably grown in recent times \cite{gazz-gru-sweers}, 
but still it is not comparable to the stage of
development of the theory concerning harmonic boundary value problems.

Scalar elliptic problems of the type
\begin{equation} \label{P.1}
 \Delta^2u = h(x,u) \quad\text{in }\Omega,
\end{equation}
where $\Omega\subset \mathbb{R}^N$, with appropriate boundary conditions 
have been studied by many authors.
When $\Omega=\mathbb{R}^3$, $h(x,u)=-u^{-q}$, with $q>0$, problem \eqref{P.1} 
is related to a fourth order analog of Yamabes equation. 
We refer to \cite{CX,IG,McKennaReichel}. In both references the authors 
studied the existence and properties of solutions.
When $\Omega=B_R$, the ball in $\mathbb{R}^N$ of radius $R$ centered at the origin,
 and $h(x,u)=\lambda\frac{f(x)}{(1-u)^q}$, problem \eqref{P.1} also 
arises in the study of MEMS. We refer to \cite{cog,GLY}, and references therein.
For general domains, and $h(x,u)=\lambda f(u)$, where the nonlinearity $f$ 
could be superlinear or singular, we refer to \cite{CEG}, where the 
regularity of the extremal solution of eigenvalue problem \eqref{P.1} is considered.

Elliptic systems of the type
\begin{equation} \label{P.2}
\begin{gathered}
 \Delta^2u + c\Delta u= f(x,u,v) \quad\text{in }\Omega,\\
  \Delta^2v + c\Delta u= g(x,u,v) \quad\text{in }\Omega,\\
\end{gathered}
\end{equation}
without singularity conditions or with appropriate singularity built into $f$ 
and $g$, $c\geq 0$ and appropriate boundary conditions have been studied 
by the some authors, see \cite{LF,JC} and references therein.
In \cite{LF} the author studied the existence result for problem \eqref{P.2},
 under Navier boundary condition,  where 
\[
f(x,u,v)=\frac{A(x)}{u^{\alpha}} + \frac{B(x)}{(u+v)^{r_1}}, \quad
g(x,u,v)=\frac{C(x)}{v^{\beta}} + \frac{D(x)}{(u+v)^{r_2}}, 
\]
 $1\leq N\leq 3$ and $r_1=r_2$. The author also proved a uniqueness result when 
$B(x)=D(x)$ and, if $A(x)=C(x)\equiv 0$, the assumption $B=D$ is not necessary 
to establish the uniqueness, see \cite[Remark 3]{LF}. The main result is 
obtained by using a version of approximating process and  Brouwer's fixed 
point theorem known as Galerking's method.
In \cite{JC} the authors used degree theory to study problem \eqref{P.2} 
with $f(x,u,v)=f(u+v)$ and $g(x,u,v)=g(u+v)$, with Dirichlet boundary 
condition and without considering any singularities.

In \cite{bdf}, the authors studied the  system
\begin{equation}\label{P.3}
\begin{gathered}
 \Delta(|\Delta u_i|^{p-2}\Delta u_i)=\lambda_iw_i(x)f_i(u_1,\ldots,u_m), 
\quad\text{in }B_1,\\
  u_i=\Delta u_i=0,\quad  x\in\partial B_1,\; i=1,\ldots,m,
\end{gathered}
\end{equation}
where $B_1$ is the unit ball in $\mathbb{R}^2$  centered at the origin, 
$p>1$ and $m\geq 1$ are integers, $w_i$ is radially symmetric,  
$f_i$ is a positive continuous function and $f_i(y_1,\ldots,y_m)$ 
may be singular at $y_i=0$. Under suitable conditions, the authors 
discuss the existence, uniqueness and dependence of solutions on the
 parameters $\lambda_i$.

In this article we generalize the result by Hernandez and Choi \cite{HC} 
for the system case and, in  cases $3<N$, the result by Faria in \cite{LF}. 
We also obtain uniqueness results in some situations that were not
 considered in \cite{LF}.

The existence of solutions for problem \eqref{P} is obtained since the 
functions $A$, $B$, $C$ and $D$ satisfy the  assumptions
\begin{equation}\label{id1}
	\begin{gathered}
 z_1=\max\big\{\min_{x \in \overline{\Omega}}\{A(x)\}, 
 \min_{x \in \overline{\Omega}}\{B(x)\}\big\}>0,\\
  z_2=\max\big\{\min_{x \in \overline{\Omega}}\{C(x)\}, 
\min_{x \in \overline{\Omega}}\{D(x)\}\big\}>0.
\end{gathered}
\end{equation}
Our main result concerning \eqref{P} is  the following.

\begin{theorem} \label{TP}
Assuming that  $\alpha,\beta,r_1,r_2\in(0,1)$ and 
$A,B,C,D\in C(\overline{\Omega})$ are nonnegative functions satisfying \eqref{id1}, 
there exists a classical solution $U=(u,v)$ of \eqref{P}.
\end{theorem}

By a (classical) solution of \eqref{P} we mean a pair of functions
 $U=(u,v) \in (C^4(\Omega) \cap C^3(\overline{\Omega}))^2$ satisfying 
the system \eqref{P}.

\begin{remark} \label{rmk1.2} \rm
In this paper we prove that if $U=(u,v)$ is a classical solution to 
problem \eqref{P}, then there exits $\delta>0$ (where $\delta$ depends 
on the sup norm $|(u, v)|_{\infty} = |u|_{\infty} + |v|_{\infty}$) so that 
$u\geq\delta\varphi_1$ and $v\geq\delta\varphi_1$, where 
$\varphi_1>0$ is the first eigenfunction of the negative Laplacian operator 
subject to zero Dirichlet boundary conditions.
\end{remark}

\begin{remark} \label{rmk1.3} \rm
 If $1\leq N\leq 3$, there exists a positive constant $\delta$ so that 
$u\geq \delta\varphi_1$ and $v\geq \delta\varphi_1$ for all classical 
solutions $(u,v)$ of \eqref{P}. Here $\delta$ does not depend on $(u,v)$, 
see Remark \ref{r3} in the Appendix.
\end{remark}

\begin{theorem}\label{TP2}
 If we assume that one of the following conditions is verified, 
then problem \eqref{P} has a unique solution.
\begin{itemize}
\item[(i)] $B\equiv 0$;
\item[(ii)] $D\equiv 0$;
 \item[(iii)] $r_1=r_2=r \in(0,1)$ and $B\equiv D$;
\item[(iv)] $r_1=r_2=r \in(0,1)$ and $A=C=0$;
\item[(v)] $r_1=r_2=r \in(0,1)$, $1\leq N\leq 3$   
and there exists a constant $\Gamma$ such that
\begin{equation*} %\label{e1}
\frac{|B(x)-D(x)|}{\varphi_1^{r+1}(x)}<\Gamma \quad\text{and}\quad
\frac{r\Gamma C^2_{\Omega}}{2\delta^{r+1}}<1, 
\end{equation*}
where $\delta$ is as in Remark \ref{r3} (see also Lemma \ref{lem1}), 
 $C_{\Omega}$ is the best constant in Sobolev embedding 
$W^{2,2}(\Omega)\cap W_0^{1,2}(\Omega)\hookrightarrow L^2(\Omega)$ (see \cite{vv}).
\end{itemize}
\end{theorem}

In the prove of Theorem \ref{TP}, one of the main results is to prove 
the existence of solutions to a family of approximate problems to problem \eqref{P}.
The proof of this result is based on the Schauder's fixed point theorem.

The organization of this article is the following. 
Section 2 contains the notation used, important lemmas that will be used, 
and the study of a family of approximate problems to problem \eqref{P}.
Section 3 is devoted to the proof of Theorem \ref{TP}. 
Section 4 is devoted to the proof of Theorem \ref{TP2}.
Section 5 is devoted to obtaining a priori estimates (in the $L^{\infty}$ sense)
on the classical solutions of \eqref{P}, in the cases $1\leq N \leq 3$.

\section{Notation and auxiliary results}

In this section we collect useful results regarding problem \eqref{P}. 
Let $(x_1,\ldots,x_m)$, $(y_1,\ldots,y_m)\in \mathbb{R}^m$. 
We use $(x_1,\ldots,x_m)\leq(y_1,\ldots,y_m)$ to denote 
$x_i\leq y_i$, $i=1,\ldots,m$.
Let $\varphi_1$ be the first eigenfunction of $(-\Delta)$ in $H_0^1(\Omega)$.  
Therefore, $\varphi_1$ satisfies
\begin{equation} \label{eigen-f}
\begin{gathered}
 -\Delta\varphi_1 =\lambda_1\varphi_1 \quad\text{in }\Omega,\\
 \varphi_1=0,  \quad\text{on } \partial\Omega,
\end{gathered}
\end{equation}
where $\lambda_1$ is the first  eigenvalue of $(-\Delta)$ with zero
 Dirichlet boundary conditions.

It is well known that $\varphi_1$ has constant sign in $\Omega$, 
so by suitable normalization we may assume $\varphi_1>0$ in $\Omega$ and  
$|\varphi_1|_{\infty}=1$. From Hopf's lemma \cite{PW}, there exists $\sigma>0$ 
such that $-\frac{\partial\varphi_1}{\partial\eta}\geq \sigma$ for all
 $x \in \partial\Omega$, where $\eta$ is the outer unit normal to
 $\partial\Omega$. Thus, $|\nabla\varphi_1(x)|\geq \sigma$ for all 
 $x \in \partial\Omega$, and there exists $c>0$ such that 
\begin{equation}\label{phi1}
c\delta_0(x)\leq \varphi_1(x)\leq\frac{1}{c}\delta_0(x),
\end{equation}
where $\delta_0(x)=\operatorname{dist}(x,\partial\Omega)$.

We denote by $G(\cdot,\cdot)$ the Green's function associated with the
 negative Laplacian operator subject to zero Dirichlet boundary conditions. 
It is known that $G$ in non-negative. If $h \in C(\overline{\Omega})$, 
the problem
 \begin{equation}\label{w}
-\Delta w = h(x)\text{ in } \Omega, \quad  w|_{\partial\Omega}=0,	
\end{equation}
has solution
 \begin{equation}\label{w2}
w(x)=\int_{\Omega}G(x,y)h(y)dy.	
\end{equation}

Now, let $\phi_0$  be the function that satisfies
\begin{equation}\label{phi0}
-\Delta\phi_0 = 1\text{ in } \Omega, \quad
 \phi_0|_{\partial\Omega}=0.	
\end{equation}
By the maximum principle we obtain $\phi_0(x)>0$ in $\Omega$.
Therefore,
\begin{gather*}
\varphi_1(x)=\lambda_1\int_{\Omega}G(x,y)\varphi_1(y)dy,\\
\phi_0(x)=\int_{\Omega}G(x,y)dy,
\end{gather*}
which, as a consequence of the normalization of $\varphi_1$, leads to
 \begin{equation}\label{w3}
\varphi_1\leq \lambda_1\phi_0.	
\end{equation}
The next lemma, due Hernandez and Choi \cite{HC}, gives an estimate 
which will be useful in proving our results.

\begin{lemma}\label{lem0}
Given $0<\xi<1$, there exists a constant $C=C(\xi)>0$, such that for all 
$x \in \Omega$,
\[
\int_{\Omega}\frac{G(x,y)}{\varphi_1(y)^{\xi}}dy \leq C(\xi).
\]
\end{lemma}

Now, for each $\epsilon \in (0,1)$ fixed, consider the auxiliary problem
\begin{equation} \label{Pepsilon}
\begin{gathered}
 \Delta^2u = \frac{A(x)}{|u|^{\alpha}+\epsilon} 
+ \frac{B(x)}{|u+v|^{r_1}+\epsilon} \quad\text{in }\Omega,\\
  \Delta^2v = \frac{C(x)}{|v|^{\beta}+\epsilon} 
+ \frac{D(x)}{|u+v|^{r_2}+\epsilon} \quad\text{in }\Omega,\\
  u=0, \quad \Delta u=0 \quad\text{on } \partial\Omega,\\
 v=0, \quad \Delta v=0 \quad\text{on } \partial\Omega,
\end{gathered}
\end{equation}
where $\Omega$ is a smooth domain in $\mathbb{R}^N$ 
($1\leq N$), $A$, $B$, $C$, $D \in C(\overline{\Omega})$ are nonnegative 
functions satisfying \eqref{id1} and $\alpha, \beta, r_1, r_2 \in (0,1)$.

The following result will be used to assure us, under some suitable 
assumptions on the nonlinearity, the uniqueness result.

\begin{lemma}\label{lem1}
Suppose that \eqref{id1} holds. Let $(u_{\epsilon}, v_{\epsilon})$,
 $\epsilon \in (0,1)$, a classical solution of \eqref{Pepsilon}. 
If there exists $K>0$, independent of $\epsilon$, such that 
$|(u_{\epsilon}, v_{\epsilon})|_{\infty} 
= |u_{\epsilon}|_{\infty} + |v_{\epsilon}|_{\infty}\leq K$, 
then there exist positive constants $\delta_1$ and $\delta_2$ 
(independent of $\epsilon$) such that
\[
(u_{\epsilon},v_{\epsilon}) \geq (\delta_1\varphi_1, \delta_2\varphi_1).
\]
\end{lemma}

\begin{proof}
Let $(u_{\epsilon}, v_{\epsilon})$ be a solution of \eqref{Pepsilon}, 
$\epsilon \in (0,1)$, $K>0$ such that 
$|(u_{\epsilon}, v_{\epsilon})|_{\infty} 
 = |u_{\epsilon}|_{\infty} + |v_{\epsilon}|_{\infty}\leq K$,
$A_0=\min_{x \in \overline{\Omega}}{A(x)}$, 
$B_0=\min_{x \in \overline{\Omega}}B(x)$, 
$C_0=\min_{x \in \overline{\Omega}}C(x)$, 
$D_0=\min_{x \in \overline{\Omega}}D(x)$. 
Let $(\omega_1, \omega_2)= (u_{\epsilon} -\delta_1\varphi_1, v_{\epsilon}
 -\delta_2\varphi_1)$, where
\begin{gather*}
0<\delta_1 < \frac{1}{\lambda_1^2}
\Big( \frac{A_0}{K^{\alpha}+1}+\frac{B_0}{(2K)^{r_1}+1}\Big), \\
0<\delta_2 < \frac{1}{\lambda_1^2}\Big( \frac{C_0}{K^{\beta}+1}
+\frac{D_0}{(2K)^{r_2}+1}\Big).
\end{gather*}]
The choice of  $\delta_1,\;\delta_2$ is always possible by \eqref{id1}. 
Then $(\Delta^2\omega_1, \Delta^2\omega_2) \geq (0,0)$ in $\Omega$, 
and $\omega_1=\omega_2=\Delta\omega_1=\Delta\omega_2=0$ on $\partial\Omega$. 
By using the Maximum Principle, we obtain $(\omega_1, \omega_2)\geq (0,0)$, and so
\[
(u_{\epsilon},v_{\epsilon}) \geq (\delta_1\varphi_1, \delta_2\varphi_1) 
\quad \text{in }  \overline{\Omega}.
\]
\end{proof}

System \eqref{Pepsilon} can be written as the system of equations
\begin{equation} \label{Pepsilon2}
\begin{gathered}
 \Delta u  + \lambda_1\,w= 0 \quad\text{in }\Omega,\\
 \Delta w + \frac{1}{\lambda_1}\frac{A(x)}{|u|^{\alpha}+\epsilon} 
+ \frac{1}{\lambda_1}\frac{B(x)}{|u+v|^{r_1}+\epsilon}=0 \quad\text{in }\Omega,\\
 \Delta v  + \lambda_1\,z= 0 \quad\text{in }\Omega,\\
  \Delta z + \frac{1}{\lambda_1}\frac{C(x)}{|v|^{\beta}+\epsilon} 
+ \frac{1}{\lambda_1}\frac{D(x)}{|u+v|^{r_2}+\epsilon}=0 \quad\text{in }\Omega,\\
  u=v=w=z=0 \quad\text{on } \partial\Omega.
\end{gathered}
\end{equation}
Let
\begin{align*}
\mathcal{A}=
\big\{& (u, w, v, z) \in (C(\overline{\Omega}))^4:
 (\tau_1\varphi_1,\tau_1\varphi_1,\tau_2\varphi_1,\tau_2\varphi_1) \\
&   \leq (u ,w , v , z) \leq(K_1, K_2,K_1, K_2)  \big\}.
\end{align*}
Let $(u,w,v,z)\in \mathcal{A}$, define
\begin{equation}\label{Tepsilon}
	T_{\epsilon}\begin{pmatrix}
u\\
w\\
v\\
z
\end{pmatrix}
=
\begin{pmatrix}
\lambda_1\int_{\Omega}G(x,y)w(y)dy\\
 \frac{1}{\lambda_1}\int_{\Omega}G(x,y)\frac{A(y)}{u(y)^{\alpha}+\epsilon}dy + \frac{1}{\lambda_1}\int_{\Omega}G(x,y)\frac{B(y)}{(u(y)+v(y))^{r_1}+\epsilon}dy\\
\lambda_1\int_{\Omega}G(x,y)z(y)dy\\
\frac{1}{\lambda_1}\int_{\Omega}G(x,y)\frac{C(y)}{v(y)^{\beta}+\epsilon}dy + \frac{1}{\lambda_1}\int_{\Omega}G(x,y)\frac{D(y)}{(u(y)+v(y))^{r_2}+\epsilon}dy
\end{pmatrix}.
\end{equation}

\begin{lemma}\label{lem2}
Suppose that \eqref{id1} holds. There exist $K_1, K_2, \tau_1$ and 
$\tau_2$ such that $T_{\epsilon}$ maps $\mathcal{A}$ into $\mathcal{A}$.
\end{lemma}

\begin{proof}
Let $A_0, B_0, C_0, D_0$ be as in Lemma \ref{lem1}, 
let $C(\alpha), C(\beta), C(r_1), C(r_2)$ be as defined in Lemma \ref{lem0},
 and define $A_{\infty}=\max_{x \in \overline{\Omega}}A(x)$, 
$B_{\infty}=\max_{x \in \overline{\Omega}}B(x)$, 
$C_{\infty}=\max_{x \in \overline{\Omega}}C(x)$, 
$D_{\infty}=\max_{x \in \overline{\Omega}}D(x)$, 
$m_0=\max_{x \in \overline{\Omega}}\phi_0(x)$.
 Choose $K_1$ such that
\begin{gather*}
\frac{A_0}{\lambda_1^2[(K_1^{\alpha})+1]} 
+ \frac{B_0}{\lambda_1^2[(2K_1)^{r_1}+1]} 
> \max\Big\{ \frac{(2A_{\infty}m_0C(\alpha))^{1/\alpha}}{K_1^{1/\alpha}}, 
\frac{(2B_{\infty}m_0C(r_1))^{1/r_1}}{K_1^{1/r_1}} \Big\},
\\
\frac{C_0}{\lambda_1^2[(K_1^{\beta})+1]} 
+ \frac{D_0}{\lambda_1^2[(2K_1)^{r_2}+1]} 
> \max\Big\{ \frac{(2C_{\infty}m_0C(\beta))^{1/\beta}}{K_1^{1/\beta}},
 \frac{(2D_{\infty}m_0C(r_2))^{1/r_2}}{K_1^{1/r_2}} \Big\}
\end{gather*}
which are always possible since $\alpha, \beta, r_1, r_2 \in (0,1)$ 
and by \eqref{id1}. Now choose $\tau_1$ and $\tau_2$ such that
\begin{align*}
&\frac{A_0}{\lambda_1^2[(K_1^{\alpha})+1]} 
+ \frac{B_0}{\lambda_1^2[(2K_1)^{r_1}+1]} \\
&> \tau_1
>\max\Big\{ \frac{(2A_{\infty}m_0C(\alpha))^{1/\alpha}}{K_1^{1/\alpha}}, 
\frac{(2B_{\infty}m_0C(r_1))^{1/r_1}}{K_1^{1/r_1}} \Big\},
\\
&\frac{C_0}{\lambda_1^2[(K_1^{\beta})+1]} 
+ \frac{D_0}{\lambda_1^2[(2K_1)^{r_2}+1]} \\
&> \tau_2
> \max\Big\{ \frac{(2C_{\infty}m_0C(\beta))^{1/\beta}}{K_1^{1/\beta}}, 
\frac{(2D_{\infty}m_0C(r_2))^{1/r_2}}{K_1^{1/r_2}} \Big\}.
\end{align*}
Then
\begin{gather*}
\frac{A_{\infty}m_0C(\alpha)}{\tau_1^{\alpha}} 
+ \frac{B_{\infty}m_0C(r_1)}{\tau_1^{r_1}}<K_1,
\\
\frac{C_{\infty}m_0C(\beta)}{\tau_2^{\beta}} 
+ \frac{D_{\infty}m_0C(r_2)}{\tau_2^{r_2}}<K_1,
\\
\tau_1 < \frac{A_0}{\lambda_1^2[(K_1^{\alpha})+1]} 
+ \frac{B_0}{\lambda_1^2[(2K_1)^{r_1}+1]}
\\
\tau_2<\frac{C_0}{\lambda_1^2[(K_1^{\beta})+1]} 
+ \frac{D_0}{\lambda_1^2[(2K_1)^{r_2}+1]}.
\end{gather*}
Finally choose 
\[
K_2=\frac{K_1}{\lambda_1m_0}.
\]
With such choices of $K_1, K_2, \tau_1$ and $\tau_2$ in $\mathcal{A}$, 
we prove that $T_{\epsilon}$ maps  $\mathcal{A}$ into $\mathcal{A}$ 
by the following calculations. Without loss of generality, we 
take $0<\epsilon < 1$.
\smallskip


\noindent\textbf{Step one.}
 Let us obtain an estimate from below for $T_{\epsilon}(u,w,v,z)$.
\begin{align*}
T_{\epsilon}\begin{pmatrix}
u\\
w\\
v\\
z
\end{pmatrix}
&=\begin{pmatrix}
\lambda_1\int_{\Omega}G(x,y)w(y)dy\\
 \frac{1}{\lambda_1}\int_{\Omega}G(x,y)\frac{A(y)}{u(y)^{\alpha}+\epsilon}dy + \frac{1}{\lambda_1}\int_{\Omega}G(x,y)\frac{B(y)}{(u(y)+v(y))^{r_1}+\epsilon}dy\\
\lambda_1\int_{\Omega}G(x,y)z(y)dy\\
\frac{1}{\lambda_1}\int_{\Omega}G(x,y)\frac{C(y)}{v(y)^{\beta}+\epsilon}dy + \frac{1}{\lambda_1}\int_{\Omega}G(x,y)\frac{D(y)}{(u(y)+v(y))^{r_2}+\epsilon}dy
\end{pmatrix}\\
&\geq  \begin{pmatrix}
\lambda_1\tau_1\int_{\Omega}G(x,y)\varphi_1(y)dy\\
 \Big(\frac{A_0}{\lambda_1[(K_1^{\alpha})+1]} 
+ \frac{B_0}{\lambda_1[(2K_1)^{r_1}+1]}\Big)\int_{\Omega}G(x,y)dy\\
\lambda_1\tau_2\int_{\Omega}G(x,y)\varphi_1(y)dy\\
\left(\frac{C_0}{\lambda_1[(K_1^{\beta})+1]} + \frac{D_0}{\lambda_1[(2K_1)^{r_2}+1]}\right)\int_{\Omega}G(x,y)dy
\end{pmatrix} \\
&\geq  \begin{pmatrix}
\lambda_1\tau_1\int_{\Omega}G(x,y)\varphi_1(y)dy\\
 \Big(\frac{A_0}{\lambda_1[(K_1^{\alpha})+1]} 
+ \frac{B_0}{\lambda_1[(2K_1)^{r_1}+1]}\Big)\phi_0(x)\\
\lambda_1\tau_2\int_{\Omega}G(x,y)\varphi_1(y)dy\\
\Big(\frac{C_0}{\lambda_1[(K_1^{\beta})+1]} 
+ \frac{D_0}{\lambda_1[(2K_1)^{r_2}+1]}\Big)\phi_0(x)
\end{pmatrix}.
\end{align*}
Using  inequality \eqref{w3}, we obtain
\[
T_{\epsilon}\begin{pmatrix}
u\\
w\\
v\\
z
\end{pmatrix}
\geq  \begin{pmatrix}
\tau_1\varphi_1\\
 \Big(\frac{A_0}{\lambda_1^2[(K_1^{\alpha})+1]} 
+ \frac{B_0}{\lambda_1^2[(2K_1)^{r_1}+1]}\Big)\varphi_1(x)\\
\tau_2\varphi_1\\
\Big(\frac{C_0}{\lambda_1^2[(K_1^{\beta})+1]} 
+ \frac{D_0}{\lambda_1^2[(2K_1)^{r_2}+1]}\Big)\varphi_1(x)
\end{pmatrix} 
\geq \begin{pmatrix}
\tau_1\varphi_1\\
\tau_1\varphi_1\\
\tau_2\varphi_1\\
\tau_2\varphi_1
\end{pmatrix}.
\]


\noindent\textbf{Step two.} 
Let us obtain an estimate from above for 
$T_{\epsilon}(u,w,v,z)$.
\begin{align*}
T_{\epsilon}\begin{pmatrix}
u\\
w\\
v\\
z
\end{pmatrix}
&= \begin{pmatrix}
\lambda_1\int_{\Omega}G(x,y)w(y)dy\\
 \frac{1}{\lambda_1}\int_{\Omega}G(x,y)\frac{A(y)}{u(y)^{\alpha}
+\epsilon}dy + \frac{1}{\lambda_1}\int_{\Omega}G(x,y)\frac{B(y)}{(u(y)+v(y))^{r_1}
+\epsilon}dy\\
\lambda_1\int_{\Omega}G(x,y)z(y)dy\\
\frac{1}{\lambda_1}\int_{\Omega}G(x,y)\frac{C(y)}{v(y)^{\beta}+\epsilon}dy
 + \frac{1}{\lambda_1}\int_{\Omega}G(x,y)\frac{D(y)}{(u(y)+v(y))^{r_2}+\epsilon}dy
\end{pmatrix} \\
&\leq  \begin{pmatrix}
\lambda_1K_2\int_{\Omega}G(x,y)dy\\
 \frac{A_{\infty}}{\lambda_1}\int_{\Omega}G(x,y)
 \frac{1}{(\tau_1\varphi_1)^{\alpha}}dy 
 + \frac{B_{\infty}}{\lambda_1}\int_{\Omega}G(x,y)\frac{1}
 {(\tau_1+\tau_2)^{r_1}\varphi^{r_1}}dy\\
\lambda_1K_2\int_{\Omega}G(x,y)dy\\
\frac{C_{\infty}}{\lambda_1}\int_{\Omega}G(x,y)\frac{1}{(\tau_2\varphi_1)^{\beta}}dy 
+ \frac{D_{\infty}}{\lambda_1}\int_{\Omega}G(x,y)
\frac{1}{(\tau_1+\tau_2)^{r_2}\varphi^{r_2}}dy
\end{pmatrix} \\
&\leq  \begin{pmatrix}
\lambda_1K_2m_0\\
 \frac{A_{\infty}}{\lambda_1\tau_1^{\alpha}}\int_{\Omega}G(x,y)
\frac{1}{\varphi_1^{\alpha}}dy + \frac{B_{\infty}}{\tau_1^{r_1}\lambda_1}
\int_{\Omega}G(x,y)\frac{1}{\varphi^{r_1}}dy\\
\lambda_1K_2m_0\\
\frac{C_{\infty}}{\lambda_1\tau_2^{\beta}}
\int_{\Omega}G(x,y)\frac{1}{\varphi_1^{\beta}}dy
 + \frac{D_{\infty}}{\lambda_1\tau_2^{r_2}}\int_{\Omega}G(x,y)
\frac{1}{\varphi^{r_2}}dy\\
\end{pmatrix}.
\end{align*}
Using  Lemma \ref{lem0}, we obtain
\[
T_{\epsilon}\begin{pmatrix}
u\\
w\\
v\\
z
\end{pmatrix}
\leq  \begin{pmatrix}
\lambda_1K_2m_0\\
 \frac{A_{\infty}C(\alpha)}{\lambda_1\tau_1^{\alpha}} 
+ \frac{B_{\infty}C(r_1)}{\lambda_1\tau_1^{r_1}}\\
\lambda_1K_2m_0\\
\frac{C_{\infty}C(\beta)}{\lambda_1\tau_2^{\beta}} 
+ \frac{D_{\infty}C(r_2)}{\lambda_1\tau_2^{r_2}}
\end{pmatrix} 
\leq  \begin{pmatrix}
\lambda_1K_2m_0\\
 \frac{K_1}{m_0\lambda_1}\\
\lambda_1K_2m_0\\
\frac{K_1}{m_0\lambda_1}
\end{pmatrix}
\leq\begin{pmatrix}
K_1\\
 K_2\\
K_1\\
K_2
\end{pmatrix}.
\]
Thus $T_{\epsilon}$ maps  $\mathcal{A}$ into $\mathcal{A}$ which 
complete the proof of Lemma \ref{lem2}.
\end{proof}

\section{Proof of Theorem \ref{TP}}

This proof will be done by means of Schauder fixed point theorem.
By Lemma \ref{lem2},  we can define $T_{\epsilon}:\mathcal{A}\to\mathcal{A}$.
Notice that $\mathcal{A}$ is closed and convex.
Now, we want to prove that the map $T_{\epsilon}$ is compact. 
In fact, let $(u\; w\;v\;z)\in \mathcal{A}$. Considering system 
\eqref{Pepsilon2}, since 
$$
\Lambda=\begin{pmatrix} w\\
 \frac{A}{|u|^{\alpha}+\epsilon} +\frac{B}{|u+v|^{r_1}+\epsilon}\\
 z\\ 
\frac{C}{|v|^{\beta}+\epsilon} + \frac{D}{|u+v|^{r_2}+\epsilon}
\end{pmatrix}
$$
belongs to $(C(\overline{\Omega}))^4$, then $\Lambda\in (L^p(\Omega))^4$ 
for any $1<p<\infty$. By using elliptic estimates \cite{Ag}, we obtain 
$T_{\epsilon}(u, w, v, z)\in (W^{2,p}(\Omega))^4$, for any $1<p<\infty$. 
The Sobolev-Morrey's imbedding theorem entails 
$T_{\epsilon}(u, w, v, z)\in (C^{1+\rho}(\Omega))^4$ for any $0<\rho<1$. 
This implies that $T_{\epsilon}$ is compact.

Now, rely on Schauder's fixed point theorem we obtain the existence of 
a fixed point 
$(u_{\epsilon}, w_{\epsilon}, v_{\epsilon},z_{\epsilon})\in (C^{1+\rho}(\Omega))^4$ 
of $T_{\epsilon}$. That is,
\begin{equation}\label{lim}
\begin{pmatrix}
u_{\epsilon}\\
w_{\epsilon}\\
v_{\epsilon}\\
z_{\epsilon}
\end{pmatrix}
=\begin{pmatrix}
\lambda_1\int_{\Omega}G(x,y)w_{\epsilon}(y)dy\\
 \frac{1}{\lambda_1}\int_{\Omega}G(x,y)
\frac{A(y)}{|u_{\epsilon}(y)|^{\alpha}+\epsilon}dy 
+ \frac{1}{\lambda_1}\int_{\Omega}G(x,y)
\frac{B(y)}{|u_{\epsilon}(y)+v_{\epsilon}(y)|^{r_1}+\epsilon}dy\\
\lambda_1\int_{\Omega}G(x,y)z_{\epsilon}(y)dy\\
\frac{1}{\lambda_1}\int_{\Omega}G(x,y)
\frac{C(y)}{|v_{\epsilon}(y)|^{\beta}+\epsilon}dy 
+ \frac{1}{\lambda_1}\int_{\Omega}G(x,y)
\frac{D(y)}{|u_{\epsilon}(y)+v_{\epsilon}(y)|^{r_2}+\epsilon}dy
\end{pmatrix},
\end{equation}
where $A,B,C,D\in C(\overline{\Omega})$. By bootstrap arguments we obtain 
$(u_{\epsilon},v_{\epsilon})\in (C^{4+\rho}(\overline{\Omega}))^2$ and 
$(w_{\epsilon},z_{\epsilon})\in (C^{2+\rho}(\overline{\Omega}))^2$.
By compactness results, we can extract convergent  subsequences in 
$C^{2+\widetilde{\rho}}(\overline{\Omega})$, namely 
$(u_n)$, $(w_n)$, $(v_n)$, $(z_n)$,  of $(u_{\epsilon})$, $(w_{\epsilon})$, 
$(v_{\epsilon})$, $(z_{\epsilon})$, respectively. 
Since $(u_n,w_n,v_n,z_n)\in \mathcal{A}$, there exist $\tau_1$ and $\tau_2$, 
independent of $n$, such that
$$
(\tau_1\varphi_1,\tau_1\varphi_1,\tau_2\varphi_1,\tau_2\varphi_1)    
\leq ( u_n, w_n, v_n,z_n). 
$$
By Lemma \ref{lem0}, we have 
\begin{equation}\label{11}
\frac{|G(x,y)|}{\varphi_{1}^{s}} \in L^{1}(\Omega),
\end{equation}
for all $s\in(0,1)$.
Using the Theorem
of the Dominated Convergence in \eqref{lim}, we obtain
\begin{equation}
\begin{pmatrix}
u\\
w\\
v\\
z
\end{pmatrix}
=
\begin{pmatrix}
\lambda_1\int_{\Omega}G(x,y)w(y)dy\\
 \frac{1}{\lambda_1}\int_{\Omega}G(x,y)\frac{A(y)}{|u(y)|^{\alpha}}dy 
+ \frac{1}{\lambda_1}\int_{\Omega}G(x,y)\frac{B(y)}{|u(y)+v(y)|^{r_1}}dy\\
\lambda_1\int_{\Omega}G(x,y)z(y)dy\\
\frac{1}{\lambda_1}\int_{\Omega}G(x,y)\frac{C(y)}{|v(y)|^{\beta}}dy 
+ \frac{1}{\lambda_1}\int_{\Omega}G(x,y)\frac{D(y)}{|u(y)+v(y)|^{r_2}}dy
\end{pmatrix}.
\end{equation}
Therefore, according to our construction we have a classical solution 
$(u,v)\in (C^2(\overline{\Omega})\cap C^4(\Omega))^2$. 
To show that $(u,v)\in (C^3(\overline{\Omega}))^2$, we can follow similar 
ideas of \cite{HC}, and this completes the proof of the existence.


\section{Proof of Theorem \ref{TP2}}

Before starting the proof of the uniqueness, let us discuss the hypothesis 
(v). Note that if $B=D\pm\varepsilon\Phi\varphi_1^{r}$ is so that
 $\Phi\in C^2_0(\overline{\Omega})$ is a solution of \eqref{sec}, 
then there exists $\varepsilon_0$ such that $|B-D|$ satisfies the hypothesis
 (v) for all $\varepsilon\in (0,\varepsilon_0)$. 
In fact, given $f\in C^{\infty}(\overline{\Omega})$, let 
$\zeta\in C^2_0(\overline{\Omega})$ be the solution of
\begin{equation}\label{sec}
\begin{gathered}
 -\Delta \zeta = f \quad\text{in }\Omega,\\
  \zeta=0,  \quad\text{on } \partial\Omega.\\
\end{gathered}
\end{equation}
By Calder\'{o}n-Zygmund estimates (see \cite{GT}),
$$
\|\zeta\|_{W^{2,p'}}\leq C\|f\|_{L^{p'}}.
$$
Since $p'>N$, it follows from Morrey's imbedding that
$$
\|\zeta/\delta_0\|_{L^{}\infty}
\leq C(\|\zeta\|_{L^\infty}+\|\nabla\zeta\|_{L^\infty})
\leq C\|\zeta\|_{W^{2,p'}}
$$ 
where $\delta_0(x)=\operatorname{dist}(x,\partial\Omega)$.
By using \eqref{phi1}, there exists $\varepsilon_0>0$ such that 
$|B-D|$ satisfies the hypothesis $(v)$ for all 
$\varepsilon\in (0,\varepsilon_0)$.


Let $\delta=\min\{\delta_1,\delta_2\}$, where 
$\delta_1,\delta_2$ are given in Lemma \ref{lem1} (see also Remark \ref{r3}).

\emph{Assume condition (i) holds.}
In this proof we adapt arguments used in \cite{HC} as follows.
 Let  $U=(u,v)$ and $\widehat{U}=(\widehat{u},\widehat{v})$ be two 
classical solutions to problem \eqref{P}. By  \cite{HC} we obtain 
that $u=\widehat{u}$. Let  $\delta>0$ be such that 
$u, v,\widehat{v}\geq\delta\varphi_1$ in $\Omega$. Define  
$z_2=v-\widehat{v}$. By the Mean Value Theorem, we arrive at
$$
\Delta^2 z_2=-\beta\frac{C(x)}{\widetilde{v}^{\beta+1}}z_2-r_2
\frac{D(x)}{\widetilde{u+v}^{r_2+1}}z_2,
$$
where $\widetilde{v},\widetilde{u+v}\geq\delta\varphi_1$ in $\Omega$. 
We then multiply the previous  equations by $z_2$, respectively, and 
integrate over a smooth domain $\Omega_1$ compactly contained in $\Omega$. 
After applying the Divergence Theorem, we obtain
\begin{align*}
&\int_{\Omega_1}\big[(\Delta z_2)^2+ \beta\frac{C(x)}{\widetilde{v}^{\beta+1}}z_2^2
+r_2\frac{D(x)}{\widetilde{u+v}^{r_2+1}}z_2^2\big]dx\\
&=\int_{\partial\Omega_1}\big[\Delta z_2\frac{\partial z_2}{\partial \nu}\big]ds
-\int_{\partial\Omega_1}\big[z_2\frac{\partial\Delta z_2}{\partial \nu}\big]ds,
\end{align*}
where $\nu$ is the unit outward normal on $\partial\Omega$.

When $\Omega_1\to\Omega$, the right-hand  side of the
equation vanishes. Since $\int_{\Omega}(\Delta z_2)^2dx<\infty$, we have
$$
\int_{\Omega}\big[\beta\frac{C(x)}{\widetilde{v}^{\beta+1}}z_2^2
+r_2\frac{D(x)}{\widetilde{u+v}^{r_2+1}}z_2^2\big]dx
$$
is well defined. Hence
$$
\int_{\Omega}\big[(\Delta z_2)^2+ \beta\frac{C(x)}{\widetilde{v}^{\beta+1}}z_2^2
+r_2\frac{D(x)}{\widetilde{u+v}^{r_2+1}}z_2^2\big]dx=0
$$
which implies that $z_2=0$ in $\Omega$.

\emph{Assume condition (ii) holds.}
 The proof is similar to (i).

\emph{Assume conditions  (iii) or (iv)  holds.}
We can follow the same idea used in \cite{LF}.

\emph{Assume condition (v) holds.}
Define $u=u_1-u_2$ and $v=v_1-v_2$, where $U_1=(u_1,v_1)$ and 
$U_2=(u_2,v_2)$ are two classical solutions for problem \eqref{P}. 
Hence $u_1,u_2 \geq \delta\varphi_1$ and $v_1,v_2 \geq \delta\varphi_1$, where 
$\delta$ does not depend  of $u_i$ and $v_i$, $i=1,2$ (see Remark \ref{r3}). 
By the Mean Value Theorem,
\begin{gather*}
 \Delta^2u = -\alpha\frac{A(x)}{\overline{u}^{\alpha+1}}u
 - r\frac{B(x)}{(\overline{u+v})^{r+1}}(u+v) \quad\text{in }\Omega,\\
  \Delta^2v = -\beta\frac{C(x)}{\overline{v}^{\beta+1}}v 
- r\frac{D(x)}{(\overline{u+v})^{r+1}}(u+v) \quad\text{in }\Omega,
\end{gather*}
where $\overline{u}\geq \delta\varphi_1$, $\overline{v}\geq \delta\varphi_1$, 
$\overline{u+v}\geq \delta\varphi_1$ in $\Omega$. We then multiply the previous 
equation by $u$ and $v$, respectively, and integrate over a smooth domain 
$\Omega_1$ compactly contained in $\Omega$. After applying the Divergence 
Theorem twice, we obtain
\begin{align*}
 &\int_{\Omega_1}\big[(\Delta u)^2 + \alpha\frac{A(x)}{\overline{u}^{\alpha+1}}u^2 
+ r\frac{B(x)}{(\overline{u+v})^{r+1}}(u^2+uv)\big] \\
 &= \int_{\partial\Omega_1}\big[\Delta u\frac{\partial\,u}{\partial\eta}\big]ds
 - \int_{\partial\Omega_1}\big[u\frac{\partial\Delta u}{\partial\eta}\big]ds,\\
 &\int_{\Omega_1}\big[(\Delta v)^2 + \beta\frac{C(x)}{\overline{v}^{\beta+1}}v^2 
+ r\frac{D(x)}{(\overline{u+v})^{r+1}}(v^2+uv)\big] \\
 &= \int_{\partial\Omega_1}\big[\Delta v\frac{\partial\,v}{\partial\eta}\big]ds 
- \int_{\partial\Omega_1}\big[v\frac{\partial\Delta v}{\partial\eta}\big]ds,
\end{align*}
where $\eta$ is the unit outward normal on $\partial\Omega$. 
By adding the last two equations and using the Holder's inequality, we obtain
\begin{align*}
 &\int_{\Omega_1}\big[(\Delta u)^2 + (\Delta v)^2 
+ \alpha\frac{A(x)}{\overline{u}^{\alpha+1}}u^2 
+ \beta\frac{C(x)}{\overline{v}^{\beta+1}}v^2\big]dx \\
&\leq  \frac{r}{2}\int_{\Omega_1}\frac{|B(x) - D(x)|}
 {(\overline{u+v})^{r+1}}(u^2+v^2)dx\\
&\quad + \int_{\partial\Omega_1}\big[\Delta u\frac{\partial\,u}{\partial\eta} 
 - u\frac{\partial\Delta u}{\partial\eta} 
 + \Delta v\frac{\partial\,v}{\partial\eta} 
 - v\frac{\partial\Delta v}{\partial\eta}\big]ds\\
&\leq  \frac{r\Gamma}{2\delta^{r+1}}\int_{\Omega_1}(u^2+v^2)dx\\
 & + \int_{\partial\Omega_1}\big[\Delta u\frac{\partial\,u}{\partial\eta}
 - u\frac{\partial\Delta u}{\partial\eta}
 + \Delta v\frac{\partial\,v}{\partial\eta}
 - v\frac{\partial\Delta v}{\partial\eta}\big]ds.
\end{align*}
Taking to the limit as $\Omega_1 \to \Omega$ the right-hand sides of the 
equations approach
\[
\frac{r\Gamma }{2\delta^{r+1}}\int_{\Omega}( u^2+ v^2)dx.
\]
Since $\int_{\Omega}(\Delta u)^2dx$,
 $\int_{\Omega}(\Delta v)^2dx$ and
$\frac{r\Gamma }{2}\int_{\Omega}( u^2+ v^2)dx < \infty$, it follows that
\[ 
\int_{\Omega}\big[\alpha\frac{A(x)}{\overline{u}^{\alpha+1}}u^2 
+ \beta\frac{C(x)}{\overline{v}^{\beta+1}}v^2\big]dx
\]
is well defined. Hence, by Sobolev embedding 
$W^{2,2}(\Omega)\cap W_0^{1,2}(\Omega)\hookrightarrow L^2(\Omega)$ 
(see \cite{vv}) we have
\[ 
\int_{\Omega}\big[(\Delta u)^2 + (\Delta v)^2 \big]dx 
\leq \frac{r\Gamma C^2_{\Omega}}{2\delta^{r+1}}
\int_{\Omega}((\Delta u)^2+(\Delta v)^2)dx.
\]
Therefore,
\[
\big(1-\frac{r\Gamma C^2_{\Omega}}{2\delta^{r+1}}\big) 
\int_{\Omega}[(\Delta u)^2 + (\Delta v)^2 ]dx \leq 0
\]
Since $(1-\frac{r\Gamma C^2_{\Omega}}{2\delta^{r+1}})>0$, we conclude 
that $u=v=0$ in $\Omega$.

\section{Appendix}

\subsection*{A priori estimates}
This section is devoted to prove a priori estimate for a classical solution 
of equation \eqref{P}, in the case $1\leq N\leq 3$.

Let  $U=(u,v) \in (C^4(\Omega) \cap C^3(\overline{\Omega}))^2$ ($u,v>0$ 
in $\Omega$) be a classical solution for system \eqref{P}. 
Multiplying the first equation of \eqref{P} by $u$ and integrating by parts in
$\Omega$, we have
\begin{equation}
 \int_{\Omega}(\Delta u)^2dx  \leq  
A_{\infty}\int_{\Omega}u^{1-\alpha}dx + B_{\infty}\int_{\Omega}u^{1-r_1}dx.
\end{equation}

Since $\frac{2}{1-\alpha}, \frac{2}{1-r_1}>1$ and $u \in L^p(\Omega)$, 
for each $p\geq 1$,  by the Young's inequality we obtain, for each 
$\epsilon>0$, that
\begin{gather*}
\int_{\Omega}u^{1-\alpha}dx 
\leq \epsilon\int_{\Omega}u^{2}dx 
+ \frac{|\Omega|}{\frac{2}{1+\alpha}(\frac{2}{1-\alpha})^{\frac{1-\alpha}{1+\alpha}}
\epsilon^{\frac{1-\alpha}{1+\alpha}}},
\\
\int_{\Omega}u^{1-r_1}dx \leq \epsilon\int_{\Omega}u^{2}dx 
+ \frac{|\Omega|}{\frac{2}{1+r_1}(\frac{2}{1-r_1})^{\frac{1-r_1}{1+r_1}}
\epsilon^{\frac{1-r_1}{1+r_1}}}.
\end{gather*}
Therefore,
\[
 \int_{\Omega}(\Delta u)^2 dx \leq   
2\epsilon\max\{A_{\infty}, B_{\infty}\}\int_{\Omega}u^{2} dx 
+C_1(\alpha,r_1,\epsilon),
\]
where 
$$
C_1(\alpha,r_1,\epsilon):= \frac{|\Omega|}{\frac{2}{1+\alpha}
(\frac{2}{1-\alpha})^{\frac{1-\alpha}{1+\alpha}}
\epsilon^{\frac{1-\alpha}{1+\alpha}}}+\frac{|\Omega|}{\frac{2}{1+r_1}
(\frac{2}{1-r_1})^{\frac{1-r_1}{1+r_1}}\epsilon^{\frac{1-r_1}{1+r_1}}}.
$$
Using the Sobolev embedding 
$W^{2,2}(\Omega)\cap W_0^{1,2}(\Omega)\hookrightarrow L^2(\Omega)$ (see \cite{vv}), 
we obtain
\[
\big(1- 2\epsilon C_{\Omega}^2\max\{A_{\infty}, B_{\infty}\}\big) 
\int_{\Omega}(\Delta u)^2 dx \leq  C_1(\alpha,r_1,\epsilon).
\]
Taking $\epsilon>0$ small enough so that 
$(1- 2\epsilon C_{\Omega}^2\max\{A_{\infty}, B_{\infty}\})>0$, 
we  obtain $\|u\|_{W^{2,2}(\Omega)} \leq C$.
 In a similar way we obtain $\|v\|_{W^{2,2}(\Omega)} \leq C$.
 Hence, since $1\leq N\leq 3$, by Sobolev embedding 
$W^{2,2}(\Omega)\hookrightarrow L^{\infty}(\Omega)$, there exists $K>0$ 
depending on $(\alpha,\beta,r_1,r_2, \Omega, A, B, C, D)$  such that
\[
\|u\|_{L^{\infty}(\Omega)} + \|v\|_{L^{\infty}(\Omega)}\leq K.
\]

\begin{remark}\label{r3} \rm
Note that if $1\leq N\leq 3$, it follows from the previous discussion
 that all classical solution $U=(u,v)$ of \eqref{P} is bounded in 
$L^{\infty}$ sense. By Lemma \ref{lem1}, there exists a positive 
constant $\delta$ (independent of $u$ and $v$) such that
$(u,v)\geq (\delta\phi_1,\delta\phi_2)$.
\end{remark}

\begin{thebibliography}{00}

\bibitem{Ag} S. Agmon;
\emph{The $L_p$ approach to the Dirichlet problem. I. Regularity theorems}, 
Ann. Scuola Norm. Sup. Pisa (3) \textbf{13} 1959, 405--448.


\bibitem{bdf} J. Barrow, R. DeYeso, L. Kong, F. Petronella;
\emph{Positive radially symmetric solution for a system of quasilinear 
biharmonic equations in the plane}.
Electron. J. Differential Equations, 2015 (2015), No. 30, 11 pp.

\bibitem{cog} D. Cassani, J. M. B. do O,   N. Ghoussoub;
\emph{On a Fourth Order Elliptic Problem with a Singular Nonlinearity}. 
Advanced Nonlinear Studies, \textbf{09} (2009), 177--197.

\bibitem{cmm} T. Chan, A. Marquina, P. Mulet;
\emph{Higher-order total variation-based image restoration}, SIAM
J. Sci. Comput. 22 (2000), 503-516.

\bibitem{CX} Y. S. Choi,  X. Xu;
\emph{Nonlinear biharmonic equations with negative exponents}, 
J. Differential Equations \textbf{246} (2009) 216--234.

\bibitem{CEG} C. Cowan, P. Esposito, N. Ghoussoub, Nassif;
\emph{Regularity of extremal solutions in fourth order nonlinear eigenvalue
 problems on general domains}, 
Discrete Contin. Dyn. Syst. 28 (2010), no. 3, 1033-1050.

\bibitem{LF} L. F. O. Faria;
\emph{Existence and uniqueness of positive solutions for singular 
biharmonic elliptic systems}, Dynamical Systems, Differential Equations 
and Applications, AIMS Proceedings, 2015, (Madrid, Spain, 2014), pp. 400--408.

\bibitem{FM} W. Fulks J. S. Maybee;
\emph{A singular non-linear equation}, 
Osaka Math. J., \textbf{12} (2010), 3295--3318.

\bibitem{gazz-gru-sweers} F. Gazzola, H. Grunau, G. Sweers;
\emph{Polyharmonic boundary value problems. Positivity preserving
and nonlinear higher order elliptic equations in bounded domains}, 
Lecture Notes in Mathematics 1991. Springer-Verlag, Berlin, Germany (2010).

\bibitem{Ghergu} M. Ghergu;
\emph{A biharmonic equation with singular nonlinearity}, 
Proc. Edinburgh Math. Soc. 55 (2012), 155-166.

\bibitem{GT} D. Gilbarg, N. S. Trudinger;
\emph{Elliptic partial differential equations of second order}, 
Grundlehren der Mathematischen Wissenschaften, vol. 224, Springer-Verlag, 
Berlin, 1983.

\bibitem{GLY} Z. Guo, B. Lai, D. Ye;
\emph{Revisiting the biharmonic equation modelling electrostatic actuation 
in lower dimensions}. Proc. Amer. Math. Soc. \textbf{142} (2014), 2027--2034.

\bibitem{HC} G. L. Hernandez, Y. Choi;
\emph{Existence of solutions in a singular biharmonic nonlinear problem}, 
Proceedings of the Edinburgh Mathematical Society \textbf{36} (1993), 537--546.

\bibitem{IG} I. Guerra;
\emph{A note on nonlinear biharmonic equations with negative exponents}, 
J. Differential Equations 253 (2012), no. 11, 3147-3157.

\bibitem{JC} T. Jung, Q. Choi;
\emph{Existence of nontrivial solutions of the nonlinear biharmonic system}, 
Korean J. Math. \textbf{16} (2008) 135--143.

\bibitem{LazerMcKenna} A. C. Lazer, P. J. McKenna;
\emph{Large amplitude periodic oscillations in suspension bridges:, 
some new connections with nonlinear analysis}. SIAM Review. 32, 537--578 (1990).

\bibitem{McKennaReichel} P.J. McKenna, W. Reichel;
\emph{Radial solutions of singular nonlinear biharmonic equations and applications
 to conformal geometry},
Electron. J. Differential Equations 37 (2003), 1-13.

\bibitem{PW} M. H. Protter, H. F. Weinberg;
\emph{Maximum Principles in Differential Equations}, 
Prentice Hall, Englewood Cliffs, NJ, 1967.

\bibitem{vv} R. C. A. M. Van der Vorst;
\emph{Best constant for the embedding of the space $H^2(\Omega)
\cap H_0^1(\Omega)$ into $L^{2N/N-4}(\Omega)$},
Differential Integral Equations, \textbf{6} (1993),  259--276.

\bibitem{yk} Y. L. You, M. Kaveh;
\emph{Fourth-order partial differential equations for noise removal}, IEEE
Trans. Image Process. 9 (2000), 1723-1730.

\end{thebibliography}

\end{document}