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

\AtBeginDocument{{\noindent\small
Variational and Topological Methods:
Theory, Applications, Numerical Simulations, and Open Problems (2012).
{\em Electronic Journal of Differential Equations},
Conference 21 (2014),  pp. 197--221.
ISSN: 1072-6691.  http://ejde.math.txstate.edu,
http://ejde.math.unt.edu \newline ftp ejde.math.txstate.edu}
\thanks{\copyright 2014 Texas State University - San Marcos.}
\vspace{9mm}}

\begin{document} \setcounter{page}{197}
\title[\hfilneg EJDE-2014/Conf/21 \hfil Solutions for nonlinear parabolic systems]
{Strong bounded solutions for nonlinear parabolic systems}

\author[N. Mavinga, M. N. Nkashama \hfil EJDE-2014/Conf/21\hfilneg]
{Nsoki Mavinga, Mubenga N. Nkashama}  % in alphabetical order

\address{Nsoki Mavinga \newline
Department of Mathematics and Statistics, Swarthmore College,
Swarthmore, PA 19081, USA}
\email{mavinga@swarthmore.edu}

\address{Mubenga N. Nkashama  \newline
 Department of Mathematics, University of Alabama at Birmingham,
Birmingham, AL 35294-1170, USA}
\email{nkashama@math.uab.edu}

\thanks{Published February 10, 2014.}
\subjclass[2000]{35K55, 35K60}
\keywords{Nonlinear parabolic systems; nonlinear boundary
conditions; \hfil\break\indent 
Carath\'eodory functions;  quasimonotone properties;
sub and supersolutions}

\begin{abstract}
  In  this article we  study the existence of  strong bounded solutions for
  nonlinear parabolic systems on a domain which is bounded in space
  and unbounded in time (namely the entire real line). We use nonlinear
  iteration arguments combined with some a priori estimates to derive
  the existence results.  We also provide conditions under which we have
  a positive solution. Some examples are given to illustrate the results.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{corollary}[theorem]{Corollary}
\newtheorem{definition}[theorem]{Definition}
\newtheorem{example}[theorem]{Example}
\allowdisplaybreaks


\section{Introduction}

We consider the nonlinear parabolic system
\begin{equation}\label{Syeq1}
\begin{gathered}
 \frac{\partial u}{\partial t}(x,t)-L_1u (x,t)= f_1(x,t,u,v)\quad
  \text{in }\Omega\times {\mathbb{R}},\\
\frac{\partial v}{\partial t}(x,t)-L_2v (x,t)= f_2(x,t,u,v)\quad
\text{in } \Omega\times {\mathbb{R}},\\
\mathcal{B}_1 u(x,t) = g_1(x,t,u,v) \quad\text{on } \partial
\Omega\times {\mathbb{R}},\\
\mathcal{B}_2 v(x,t) = g_2(x,t,u,v)\quad\text{on } \partial
\Omega\times {\mathbb{R}},\\
 \sup_{\Omega\times
{\mathbb{R}}}\{|u(x,t)|, |v(x,t)|\}<\infty,
\end{gathered}
\end{equation}
where $\Omega$ is a bounded, open  and connected subset of
${\mathbb{R}^N}$ with smooth boundary $\partial \Omega$ and closure
$\overline{\Omega}$. We suppose that $L_k$  are  second
order,  uniformly elliptic differential operators with time-dependent
coefficients and $\mathcal{B}_{k}$ are  linear first-order boundary
operators which are either Dirichlet, Neumann or regular oblique type.
We suppose that the coefficients of the operators $L_k$ and $\mathcal{B}_k$ are
measurable and bounded. The reaction and the
boundary nonlinearities $f_k$ and $g_k$  are, say,
Carath\'eodory functions.

 We are interested in bounded solutions existing for all time. Steady-state,
time-periodic solutions and (bounded) attractors as well as almost-periodic
solutions are only a few examples of solutions existing for all times,
see e.g. \cite{CC09,DO00,FM90,MP09,SZ94}.
The study of nonlinear parabolic systems for large time have important
applications in ecology. Full bounded solutions are thus important in both
backward and forward dynamics.
Many papers have been devoted to the study of nonlinear parabolic systems
with given initial conditions. For some recent results in this direction and
bibliography we refer to \cite{HA76S,HA88, AL82S,P92S,CC03S, FM90, MP09, SZ94},
and others. However, to the best of  our knowledge, not much seems to be done
for system \eqref{Syeq1} for the case of nonlinear boundary conditions and
in the unbounded time-domain (namely, the entire
real line). A few results were obtained in the scalar case by the
authors \cite{NM11,NM10}. Since we are dealing with nonlinear parabolic
systems with nonlinear boundary conditions and without initial conditions,
many of the tools used for compact or semi-infinite time interval such as the
maximum principle and the fixed-point results are not directly applicable.
Thus, the need to develop new tools for studying the problem. Moreover,
we deduce the comparison principle which is valid on the entire real line
in time and used some nonlinear iteration arguments to obtain the
existence results.

  The paper is organized as follows.  In Section 2, we formulate general
assumptions which are needed throughout, and state our main results
concerning the existence of bounded solutions
existing for all times for nonlinear systems with (possibly)
nonlinear boundary conditions. We assume
that the nonlinearities  in the reaction and on the boundary satisfy
some growth conditions.
In Section 3, we state  some results on (scalar) linear parabolic
equations which are needed in the proof of our main results.
In Section 4, we prove the main results.  We conclude the paper with
some examples which illustrate our results.

\section{Assumptions and Main results}

All functions in this paper will take values in ${\mathbb{R}}$
and all vector spaces are over the reals. We assume that $\Omega$ is
a bounded domain in ${\mathbb{R}^N}$ with boundary $\partial
\Omega$ and closure $\overline{\Omega}$. We assume that $\partial \Omega$
belongs to $C^{2}$; $\mu\in(0,1)$ and $p=(N+2)/(1-\mu)$.
We consider the second
order parabolic operators in ${\Omega}\times {\mathbb{R}}$ given by
\begin{equation}\label{operator1}
\frac{\partial u}{\partial t} - L_ku,
\end{equation}
where
\[
L_ku :=\sum_{i,j=1}^N a_{ij}^{(k)}(x,t)\frac{\partial^{2}
u }{\partial x_i
\partial x_j}+ \sum_{i=1}^N b_i^{(k)}(x,t)\frac{\partial u}{\partial
x_i}+c^{(k)}(x,t)u ,\quad k=1,2
\]
with symmetric positive definite
coefficient-matrices $(a_{ij}^{(k)})$. We assume that
\begin{itemize}
\item[(i)] $a_{ij}^{(k)}\in C(\overline{\Omega}\times {\mathbb{R}})\cap
L^{\infty}(\Omega\times {\mathbb{R}})$, $b_i^{(k)}, c^{(k)}\in
L^{\infty}(\Omega\times {\mathbb{R}})$.

\item[(ii)] There are constants $c_0\ge0$ and $\gamma_0>0$ such that for
a.e. $(x,t)\in \overline\Omega\times {\mathbb{R}}$, $c^{(k)}(x,t)\leq
-c_0$ and $ \sum_{i,j=1}^N a_{ij}^{(k)}(x,t)\xi_i\xi_j\geq
\gamma_0 |\xi|^2 $ for all $\xi \in  {\mathbb{R}}^N$.
\end{itemize}

Let $\epsilon$ denote a variable which assumes the values 0 and 1
only. We define the boundary operators $\mathcal{B}_{k,\epsilon}$ by
\begin{equation}\label{form2}
\mathcal{B}_{k,\epsilon}u:= \epsilon\frac{\partial u}{\partial \nu} +
\alpha^{(k)}(x,t)u,\quad k=1,2
\end{equation}
where $\alpha^{(k)}\in W_{p,\rm{loc}}^{1,1/2}({\partial\Omega}\times
{\mathbb{R}})\cap L^{\infty}(\Omega\times {\mathbb{R}})$, and for
all $(x,t)\in {\partial\Omega}\times {\mathbb{R}}$,
$\alpha^{(k)}(x,t)\geq\alpha_0 \geq 0$. The constant $\alpha_0$ is such
that $\alpha_0>0 \; \text{if}\; \epsilon=0, \; \text{and}\;
\alpha_0\geq 0 \;\text {if}\; \epsilon=1$. Moreover,
\begin{equation}\label{positivity}
c_0+ \alpha_0>0, \end{equation}
which implies that the coefficients
$c^{(k)}(x,t)$ and $\alpha^{(k)}(x,t)$ do not vanish simultaneously. Thus, for
$\epsilon=0$, $\mathcal{B}_{k,0} u$ is a Dirichlet boundary condition
whereas for $\epsilon=1$, $\mathcal{B}_{k,1}u$ corresponds to a Neumann
or a regular oblique derivative boundary condition.

In what follows,  the inequality $(u_1,v_1)\leq (u_2,v_2)$
means that $u_1\leq u_2$ and $v_1\leq v_2$. The functions $f_k, g_k$
depend in general on $u$ and $ v;$ except for the Dirichlet boundary
condition where $g_k$ are independent of $u$ and $v$; that is,
$g_{k}(x,t,u,v)=g_k(x,t)$.
The reaction  functions $f_k\in L_{\text{car}}^\infty(\Omega\times
\mathbb{R}\times\mathbb{R}\times\mathbb{R})$ are
$L^\infty$-Carath\'eodory
functions; that is,
\begin{itemize}
\item[(A1)] $f_k(\cdot,\cdot,u,v)$ is measurable; $f_k(x,t,\cdot,\cdot)$
is continuous for a.e. $(x,t)\in\Omega\times\mathbb{R}$;
and for every $r>0$ there is a function $M_{k,r}\in
L^\infty(\Omega\times\mathbb{R})$ such that
$|f_k(x,t,u,v)|\leq M_{k,r}(x,t)$ for a.e.
$(x,t)\in\Omega\times\mathbb{R}$ and all
$(u, v)\in[-r,r]\times[-r,r]$.
\end{itemize}

The boundary functions $g_k$ are also $L^\infty$-Carath\'eodory
functions; i.e, they satisfy (A1), but in addition if $\epsilon=0$
 (i.e Dirichlet boundary condition) then
$g_{k} \in W^{2-1/p,(2-1/p)/2}_{p,\rm{loc}}(\partial\Omega\times\mathbb{R})$
and if $\epsilon=1$ (i.e Neumann boundary condition) then $g_{k}(x,t, u,v)$
satisfy the Lipschitz condition in $J_1\times J_2\subset{\mathbb{R}}^2$,
uniformly in $(x,t)\in \overline{\Omega}\times {\mathbb{R}}$, where
$J_1$ and $J_2$ are closed intervals in $\mathbb{R}$; that is
\begin{itemize}
\item[(A2)] for every $J_1\times J_2\subset {\mathbb{R}}\times  {\mathbb{R}}$,
there is a constant
$\varrho_k=\varrho_k(\partial\Omega\times\mathbb{R}\times J_1\times J_2)>0$ such
that
\begin{align*}
&|g_k(x,t,u_1,v_1)-g_k(y,s,u_2, v_2)|\\
&\leq\varrho_k\left[|x-y|^2+|t-s|+|u_1-u_2|^2+|v_1-v_2|^2\right]^{1/2}
\end{align*}
for all
$(x,t,u_1,v_1), (y,s,u_2, v_2)\in\partial\Omega\times\mathbb{R}
\times J_1\times J_2$.

\item[(A3)] The vector functions $ \mathbf{f}=(f_{1},f_{2}), \mathbf{g}
=(g_{1},g_{2})$ are  quasimonotone in $J_1\times J_2$, that
is, they satisfy one of the following quasimonotonicity properties:
$(f_1, f_2),\, (g_1, g_2)$  are
quasimonotone nondecreasing (quasimonotone nonincreasing) in
$J_1\times J_2;$ i.e.,
 $$
\text{for fixed }u\in J_1,  f_1, g_1
\text{ are nondecreasing (nonincreasing) in } v\in  J_2,
$$
and
$$
\text{for fixed }v\in J_2,  f_2,
 g_2\;\text{ are  nondecreasing (nonincreasing) in } u\in J_1.
$$
\end{itemize}

We wish to emphasize the fact that this `additional' local Lipschitz
condition (A2) on the boundary nonlinearities $g_k(x,t,u,v)$ is
needed to obtain \emph{a priori} estimates for the \emph{boundary
traces} of solutions.
It ensures that the boundary
superposition (Nemytsk\v{\i}i) operator associated with the function
$g_k(x,t,\cdot,\cdot)$ maps
$W^{1-1/p,(1-1/p)/2}_{p,\rm{loc}}(\partial\Omega\times\mathbb{R})
\times W^{1-1/p,(1-1/p)/2}_{p,\rm{loc}}(\partial\Omega\times\mathbb{R})$
into
$W^{1-1/p,(1-1/p)/2}_{p,\rm{loc}}(\partial\Omega\times\mathbb{R})$;
the latter is the condition needed on the boundary data to get
strong solutions, i.e., solutions in
$W^{2,1}_{p,\rm{loc}}(\Omega\times\mathbb{R})
\times W^{2,1}_{p,\rm{loc}}(\Omega\times\mathbb{R})$. In particular, for
$g_k(x,t,u,v)=g_k(x,t)$ independent of $(u,v)$, it implies that
$g_k\in W^{1-1/p,(1-1/p)/2}_{p,\rm{loc}}(\partial\Omega\times\mathbb{R})$.

Based on the type of quasimonotonicity property, we will use the
following definitions for strong sub and supersolutions.

\begin{definition}\label{defsubsol} \rm
 A pair of functions $(\underline{u},
\underline{v})$ and $(\overline{u}, \overline{v})$
in $W^{2,1}_{p,\rm{loc}}(\Omega\times
{\mathbb{R}})\times W^{2,1}_{p,\rm{loc}}(\Omega\times
{\mathbb{R}})$ are ordered \textit{subsolution} and
\textit{supersolution} of the system \eqref{Syeq1} if
\begin{itemize}
\item[(1)] $(\underline{u}, \underline{v})\leq (\overline{u},
\overline{v})$, that is, $\underline{u}\leq \overline{u}$ and
$\underline{v}\leq \overline{v}$, and one of the following
conditions holds.
 \item[(2)]  When $(f_1,f_2) $ and $(g_1, g_2)$ are
quasimonotone nondecreasing
\begin{gather*}
\frac{\partial \overline{u}}{\partial
t}-L_1\overline{u}-f_1(x,t,\overline{u},\overline{v})\geq 0\geq
\frac{\partial \underline{u}}{\partial t}-L_1\underline{u}-
f_1(x,t,\underline{u},\underline{v})\quad\text{in } \Omega\times
{\mathbb{R}},
\\
\frac{\partial \overline{v}}{\partial
t}-L_2\overline{v}-f_2(x,t,\overline{u},\overline{v})\geq 0\geq
\frac{\partial \underline{v}}{\partial t}-L_2\underline{v}-
f_2(x,t,\underline{u},\underline{v})\quad\text{in } \Omega\times
{\mathbb{R}},
\\
\mathcal{B}_1\overline{u}- g_1(x,t,\overline{u}, \overline{v})\geq
0\geq\mathcal{B}_1\underline{u}- g_1(x,t,\underline{u},
\underline{v}) \quad\text{on } \partial\Omega\times {\mathbb{R}},
\\
\mathcal{B}_2\overline{v}- g_2(x,t,\overline{u}, \overline{v})\geq
0\geq\mathcal{B}_2\underline{v}- g_2(x,t,\underline{u},
\underline{v}) \quad\text{on }\partial \Omega\times {\mathbb{R}},
\\
\sup_{\Omega\times {\mathbb{R}}}(|\underline{u}|,|\underline{v}|,|\overline{u}|,
|\overline{v}|)<\infty.
\end{gather*}

\item[(3)] When $(f_1,f_2) $ and $(g_1, g_2)$ are
quasimonotone nonincreasing
\begin{gather*}
\frac{\partial \overline{u}}{\partial
t}-L_1\overline{u}-f_1(x,t,\overline{u},\underline{v})\geq 0\geq
\frac{\partial \underline{u}}{\partial t}-L_1\underline{u}-
f_1(x,t,\underline{u},\overline{v})\quad\text{in } \Omega\times
{\mathbb{R}},
\\
\frac{\partial \overline{v}}{\partial
t}-L_2\overline{v}-f_2(x,t,\underline{u},\overline{v})\geq 0\geq
\frac{\partial \underline{v}}{\partial t}-L_2\underline{v}-
f_2(x,t,\overline{u},\underline{v})\quad \text{in } \Omega\times
{\mathbb{R}},
\\
\mathcal{B}_1\overline{u}- g_1(x,t,\overline{u},
\underline{v})\geq 0\geq\mathcal{B}_1\underline{u}-
g_1(x,t,\underline{u}, \overline{v}) \quad\text{on } \partial\Omega\times
{\mathbb{R}},
\\
\mathcal{B}_2\overline{v}- g_2(x,t,\underline{u},
\overline{v})\geq 0\geq\mathcal{B}_2\underline{v}-
g_2(x,t,\overline{u}, \underline{v}) \quad\text{on } \partial\Omega\times
{\mathbb{R}},
\\
\sup_{\Omega\times
{\mathbb{R}}}(|\underline{u}|,|\underline{v}|,|\overline{u}|,|\overline{v}|
)<\infty.
\end{gather*}
\end{itemize}
For the rest of this article, we assume that the interval
$[\underline{u}, \overline{u}] \times [\underline{v}, \overline{v}]
\subseteq J_1\times  J_2$.
\end{definition}
Our main result for the system \eqref{Syeq1} is given by the following
theorem.

\begin{theorem}\label{Syt1}
 Assume {\rm (A1)---(A3)} are satisfied and suppose that \eqref{Syeq1} has an ordered subsolution
$(\underline{u}, \underline{v})$ and supersolution $(\overline{u},
\overline{v})$  and $(f_1, f_2)$ and $(g_1,
 g_2)$ are quasimonotone nondecreasing (quasimonotone
 nonincreasing) in $[(\underline{u}, \underline{v}), (\overline{u}, \overline{v})]$.
Then the system { \eqref{Syeq1}} has at least one solution
$(u,v) \in W^{2,1}_{p,\rm{loc}}(\Omega\times
{\mathbb{R}})\times W^{2,1}_{p,\rm{loc}}(\Omega\times
{\mathbb{R}})$ such that
$$(\underline{u},
\underline{v})\leq (u,v)\leq (\overline{u},
\overline{v})\;\text{in } \overline{\Omega}\times {\mathbb{R}}.$$
\end{theorem}

As an immediate consequence of Theorem \ref{Syt1}, we have the following
corollary on the existence of positive full bounded solutions.

\begin{corollary}[Positive Solutions] \label{c1}
Assume that  the assumptions in Theorem \ref{Syt1} are satisfied.
Suppose that either $(f_1,f_2)$ and $(g_1,  g_2)$ are quasimonotone 
nondecreasing
and
$$
f_1(x,t, 0,0)\geq 0, \quad f_2(x,t, 0, 0)\geq 0, \quad 
g_1(x,t, 0, 0)\geq 0,\quad g_2(x,t,0,0)\geq  0.
$$
When either $(f_1, f_2)$ and $(g_1,  g_2)$ are quasimonotone nonincreasing and
$$
f_1(x,t, 0, v)\geq 0, \;f_2(x,t, u, 0)\geq 0,\; g_1(x,t, 0, v)\geq 0,
 \text{ and } g_2(x,t,u,0)\geq  0.
$$
Furthermore, assume that there exists a  nonegative supersolution
$(\overline{u}, \overline{v})$ of \eqref{Syeq1}. Then the
system \eqref{Syeq1}  has a nonnegative solution
$(u,v) \in W^{2,1}_{p,\rm{loc}}(\Omega\times
{\mathbb{R}})\times W^{2,1}_{p,\rm{loc}}(\Omega\times
{\mathbb{R}})$ such that
$(u,v)\leq (\overline{u}, \overline{v})$.
\end{corollary}

Indeed, observe that $(0,0)$ is a subsolution of the
system \eqref{Syeq1}. Therefore by Theorem \ref{Syt1}, the system
\eqref{Syeq1} has a nonnegative solution.

\section{Auxiliary  results and proof of the main result}

To prove the main result stated in the previous section, we
need some auxiliary  results on scalar linear parabolic equations in
$\Omega\times {\mathbb{R}}$. We refer to \cite{NM11, NM10} for
the proof of these results.

Consider the linear boundary value problem
\begin{equation}\label{SyL1}
\begin{gathered}
 \frac{\partial u}{\partial t}-Lu= f \quad \text{in } \Omega\times {\mathbb{R}},\\
\mathcal{B}_{\epsilon}u=\varphi \quad \text{on }
\partial \Omega \times {\mathbb{R}},\\
 \sup_{\Omega\times {\mathbb{R}}}|u(x,t)|<\infty,
\end{gathered}
\end{equation}
where $L$ is  a second order, time dependent, uniformly elliptic
differential operators and $\mathcal{B}_{\epsilon}$ is a first order
boundary operator  as defined in the previous section.

\begin{proposition}[A priori estimates] \label{Syprop1}
Let $u\in W^{2,1}_{p,\rm{loc}}(\Omega\times {\mathbb{R}})$ be
(uniformly)\\ bounded at $-\infty$. Then there exists a constant $K$
such that
\begin{equation}\label{form3}
\sup_{\Omega\times
{\mathbb{R}}}|u_{\pm}|\leq K
\Big(\sup_{\Omega\times {\mathbb{R}}}\big|\big( \frac{\partial u}{\partial t}
-Lu\big)_{\pm}\big|+
\sup_{\partial\Omega\times
{\mathbb{R}}}|(\mathcal{B}_{\epsilon}u)_{\pm}|\Big);
\end{equation}
which implies that
$$
\sup_{\Omega\times {\mathbb{R}}}|u|\leq K
\Big(\sup_{\Omega\times {\mathbb{R}}}|
 \frac{\partial u}{\partial t}-Lu|+\sup_{\partial\Omega\times
{\mathbb{R}}}|\mathcal{B}_{\epsilon}u|\Big).
$$
The constant $K$ depends only on the dimension $N$, the parabolicity
constant $\gamma_0$, $\operatorname{diam}(\Omega)$, and the
$L^{\infty}$-bounds of the coefficients of the operators $L$ and
$\mathcal{B}_{\epsilon}$.
\end{proposition}

We deduce from the above proposition the following (weak) maximum
type-comparison principle.

\begin{corollary}\label{Sycor2} {\rm (Weak Maximum/Comparison Principle)}
Suppose that the conditions  of Proposition \ref{Syprop1} are met. Assume
that $\frac{\partial u}{\partial t}-Lu \geq 0$ a.e. in
$\Omega\times {\mathbb{R}}$ and that $\mathcal{B}_\epsilon u\geq 0$
on $\partial\Omega\times {\mathbb{R}}$. Then $u\geq 0$ in
${\overline{\Omega}}\times {\mathbb{R}}$.
\end{corollary}

The next proposition  deals with  the existence result for linear
parabolic equations.

\begin{proposition}\label{Sylin}
Suppose that  $f\in L^{\infty}(\Omega\times {\mathbb{R}})$ and
$\varphi \in \mathrm{W}^{2-\epsilon-\frac{1}{p},(2-\epsilon
-\frac{1}{p})/2}_{p,\rm{loc}}(\partial\Omega\times
{\mathbb{R}})\cap L^{\infty}(\partial\Omega \times {\mathbb{R}})$
with $p=\frac{N+2}{1-\mu}$.
 Then the problem  \eqref{SyL1} has a unique solution
$u\in W^{2,1}_{p,\rm{loc}}({\Omega}\times {\mathbb{R}})\cap
L^{\infty}(\Omega\times {\mathbb{R}})$.
\end{proposition}



\begin{lemma}[Interpolation inequalities] \label{Syinterp}
 Let $\Omega\times I\subset{\mathbb{R}^n}\times {\mathbb{R}}$ and
$1\leq p<\infty$. There is a constant $C>0$ such that for all
$u\in W_p^{2,1}({\Omega}\times I)$ one has
\begin{equation*}\label{Syintp1}
|u|_{W_p^{1,1/2}({\Omega}\times I)}\leq
C |u|_{W_p^{2,1}({\Omega}\times I)}^{1/2}|u|_{L^p({\Omega}\times
I)}^{1/2}.
\end{equation*}
Moreover, for every $\varepsilon>0$,
\begin{equation} \label{Syintp2}
|u|_{W_p^{1,1/2}({\Omega}\times I)}\leq
C\Big(\varepsilon |u|_{W_p^{2,1}({\Omega}\times
I)}+\frac{1}{4\varepsilon}|u|_{L^p({\Omega}\times I)}\Big).
\end{equation}
\end{lemma}

\begin{proposition}\label{uniqueness}
 Consider the nonlinear parabolic boundary value problem
 \begin{equation}\label{eq1S}
\begin{gathered}
\frac{\partial u}{\partial t}(x,t)-Lu (x,t)= f(x,t,u) \quad \text{a.e. in }
 \Omega\times{\mathbb{R}},\\
\mathcal{B}u = \varphi(x,t,u) \quad \text{a.e. on }
 \partial \Omega\times {\mathbb{R}},\\
 \sup_{\Omega\times {\mathbb{R}}}|u(x,t)|<\infty.
\end{gathered}
\end{equation}
Suppose that $\underline{u}, \overline{u}$ are ordered sub-solution and
super-solution of  \eqref{eq1S} and suppose that $f$ and
$\varphi$  are nonincreasing in $u$, for
$u\in[\underline{u},\overline{u}]$. Then  \eqref{eq1S}  has at
most one solution $u$ such that $\underline{u}\leq
u\leq\overline{u}$.
\end{proposition}


\begin{lemma} \label{transfo1}
Let $f$ satisfy {\rm (A1)}; that is, $f\in
L^\infty_{\rm{car}}(\Omega\times {\mathbb{R}}\times\mathbb{R})$.
Then, for every $r\in\mathbb{R}$ with $r>0$, there is a continuous
function $m:[-r,r]\times[-r,r]\to\mathbb{R}$ such that $m(\cdot,v)$
is nondecreasing on $[-r,r]$, $m(u,\cdot)$ is nonincreasing on
$[-r,r]$, $m(u,v)=-m(v,u)$ on $[-r,r]\times[-r,r]$, and
\begin{equation}\label{transfo2b}
\sup_{\Omega\times\mathbb{R}}|f(x,t,u)-f(x,t,v)|\le m(u,v)
\end{equation}
for all $u,v\in[-r,r]$ with $u\ge v$.
\end{lemma}


\begin{proposition}\label{Syp11}
 Let {\rm (A1)--(A3)} and the following condition hold,
\begin{itemize}
\item[(LL)] The functions $f_k$ $(k=1,2)$  satisfy the one-sided Lipschitz
condition in $J_1\times J_2\subset{\mathbb{R}}^2$,
uniformly a.e. in $(x,t)\in \overline{\Omega}\times {\mathbb{R}}$; that
is,  there are constants $\theta_k\geq 0$ such that for every
$(u_1,v), (u_2,v), (u,v_1), (u,v_2)$ in $J_1\times J_2$,
\begin{gather*}
f_1(x,t,u_1,v)-f_1(x,t,u_2,v)\geq -\theta_1(u_1-u_2),\quad \text{for }
  u_1\geq u_2 ,\\
f_2(x,t,u,v_1)-f_2(x,t,u,v_2)\geq -\theta_2(v_1-v_2), \text {for}\;
  v_1\geq v_2  .
\end{gather*}
\end{itemize}
  Suppose that \eqref{Syeq1} has an ordered subsolution
$(\underline{u}, \underline{v})$ and supersolution $(\overline{u},
\overline{v})$  and $(f_1, f_2)$ and $(g_1, g_2)$ are quasimonotone
nondecreasing (quasimonotone  nonincreasing) in
$[(\underline{u}, \underline{v}), (\overline{u}, \overline{v})]$.
Then the system  \eqref{Syeq1} has at least one solution
$(u,v)\, \in W_{p,\rm{loc}}^{2,1}(\overline{\Omega}
\times {\mathbb{R}})\times W_{p,\rm{loc}}^{2,1}
(\overline{\Omega}\times {\mathbb{R}})$ such that
$$
(\underline{u}, \underline{v})\leq (u,v)\leq (\overline{u},
\overline{v})\;\text{in } \overline{\Omega}\times {\mathbb{R}}.
$$
\end{proposition}

\begin{proof}
Let $\delta=\max\{\theta_1, \theta_2, \rho_1,\rho_2\}$.
Consider the following modified problem.
\begin{equation}\label{Syeq2}
\begin{gathered}
\frac{\partial u}{\partial t}-L_1u +\delta u= f_1(x,t,u,v)+\delta u
\quad\text{in } \Omega\times {\mathbb{R}},\\
\frac{\partial v}{\partial t}-L_2v +\delta v= f_2(x,t,u,v)+ \delta v
\quad\text{in } \Omega\times {\mathbb{R}},\\
\mathcal{B}_1 u+\delta u = g_1(x,t,u,v)+\delta u
\quad\text{on } \partial \Omega\times{\mathbb{R}},\\
\mathcal{B}_2 v +\delta v = g_2(x,t,u,v)+\delta v
\quad\text{on } \partial \Omega\times{\mathbb{R}},\\
 \sup_{\Omega\times {\mathbb{R}}}\{|u(x,t)|, |v(x,t)|\}<\infty.
\end{gathered}
\end{equation}
To prove the existence of solutions for problem \eqref{Syeq1},
it suffices  to show that the modified problem \eqref{Syeq2} has a
solution.

First, we construct a sequence $(u_n,v_n)$ from the
(linear) iteration process
\begin{equation}\label{Syeq3}
\begin{gathered}
\frac{\partial {u}_n}{\partial t}-L_1{u}_n +\delta{u}_n
= f_1(x,t,{u}_{n-1},{v}_{n-1})+\delta{u}_{n-1}\quad\text{in } \Omega\times
{\mathbb{R}},\\
\frac{\partial {v}_n}{\partial t}-L_2{v}_n +\delta{v}_n=
f_2(x,t,{u}_{n-1},{v}_{n-1})+\delta{v}_{n-1}\quad\text{in } \Omega\times
{\mathbb{R}},\\
\mathcal{B}_1 {u}_n+\delta{u}_n = g_1(x,t,{u}_{n-1},{v}_{n-1})+\delta{u}_{n-1}
\quad\text{on }
\partial \Omega\times
{\mathbb{R}},\\
\mathcal{B}_2 {v}_n +\delta{v}_n=
g_2(x,t,{u}_{n-1},{v}_{n-1})+\delta{v}_{n-1}\quad\text{on }
\partial \Omega\times
{\mathbb{R}},\\
 \sup_{\Omega\times
{\mathbb{R}}}\{|{u}_n(x,t)|, |{v}_n(x,t)|\}<\infty,
\end{gathered}
\end{equation}
where the initial iteration $(u_0,v_0)$ is determined by the
quasimonotonicity property considered.

(i) If  $f_k$ and $g_k$ are quasimonotone nondecreasing then  we
take either  $(u_0,v_0)= (\underline{u},\underline{v})$ or
 $(u_0,v_0)= (\overline{u},\overline{v})$, and we denote the sequences
 constructed from  the two initial  iterations by $ (\underline{u}_n, \underline{v}_n)$ and
$(\overline{u}_n, \overline{v}_n)$, respectively.

(ii) If $f_k$ and $g_k$ are quasimonotone nonincreasing then  we
choose either  $(u_0,v_0)= (\underline{u},\overline{v})$ or
 $(u_0,v_0)= (\overline{u},\underline{v})$, and we  denote the
 sequences constructed from the two initial  iterations by $(\underline{u}_n, \overline{v}_n)$ and
$(\overline{u}_n, \underline{v}_n)$, respectively.
For sake of discussion,  we will present the rest
of the proof for the case of quasimonotone nondecreasing functions.
Similar arguments can be used for the quasimonotone
nonincreasing case.

Observe that for each $n\in {\mathbb{N}}$, the above system consists
of two linear uncoupled problems and
$f_k(.,., u_0,v_0)\in L^{\infty}(\Omega\times \mathbb{R})$ and
$g_k(.,.,u_0,v_0) \in \mathrm{W}^{2-\epsilon-\frac{1}{p},
(2-\epsilon-\frac{1}{p})/2}_{p,\rm{loc}}\\ (\partial\Omega\times
{\mathbb{R}})\cap L^{\infty}(\partial\Omega \times {\mathbb{R}})$
whenever $(u_0,v_0)\in [\underline{u},\overline{u}] \times
[\underline{u},\overline{u}]$.  Set $(u_0,v_0)=(\underline{u},\underline{v})$
then it follows from proposition
\ref{Sylin} that  problem \eqref{Syeq3} has a unique solution
$(\underline{u}_1, \underline{v}_1)\in \mathrm{W}^{2,1}_{p,\rm{loc}}(\Omega\times
{\mathbb{R}})\cap L^{\infty}(\Omega \times {\mathbb{R}})\times
\mathrm{W}^{2,1}_{p,\rm{loc}}(\Omega\times
{\mathbb{R}})\cap L^{\infty}(\Omega \times {\mathbb{R}})$.
Moreover if $(u_0,v_0)=(\overline{u},\overline{v})$ then a similar argument
shows that problem \eqref{Syeq3} has a unique solution
$(\overline{u}_1, \overline{v}_1)\in \mathrm{W}^{2,1}_{p,\rm{loc}}(\Omega\times
{\mathbb{R}})\cap L^{\infty}(\Omega \times {\mathbb{R}})\times
\mathrm{W}^{2,1}_{p,\rm{loc}}(\Omega\times
{\mathbb{R}})\cap L^{\infty}(\Omega \times {\mathbb{R}})$.
 Furthermore, $(\underline{u},\underline{v})\leq
(\underline{u}_1, \underline{v}_1)\leq(\overline{u}_1,
\overline{v}_1)\leq(\overline{u}, \overline{v})$.

 Indeed, let $w_1=\overline{u}-\overline{u}_1$  and
$w_2=\overline{v}-\overline{v}_1$. By (A3), (LL) and the
definition of super-solution one gets
\begin{gather*}
 \frac{\partial {w_1}}{\partial t}-L_1w_1 +\delta w_1\geq 0
\quad \text{in } \Omega\times {\mathbb{R}},\\ \mathcal{B}_1 w_1+\delta w_1
\geq 0 \quad  \text{on }
\partial \Omega\times
{\mathbb{R}},\\
\sup_{\Omega\times {\mathbb{R}}}\{|w_1(x,t)|\}<\infty.
\end{gather*}
and
\begin{gather*}
\frac{\partial w_2}{\partial t}-L_2w_2
+ \delta w_2\geq 0 \quad \text{in } \Omega\times {\mathbb{R}},\\
\mathcal{B}_2 w_2 +\delta w_2 \geq 0 \quad \text{on }\partial \Omega\times
{\mathbb{R}},\\
 \sup_{\Omega\times {\mathbb{R}}}\{|w_2(x,t)|\}<\infty.
\end{gather*}
By Corollary \ref{Sycor2} it follows that $w_1\geq 0$ and
$w_2\geq 0$; that is, $\overline{u}_1\leq \overline{u}$  and
$\overline{v}_1\leq \overline{v}$.
Using the definition of subsolution we show in a similar way that
$\underline{u}\leq\underline{u}_1$ and $\underline{v}\leq
\underline{v}_1$.

We now need to prove that   $\underline{u}_1\leq\overline{u}_1$ and
$\underline{v}_1\leq\overline{v}_1$. Let
$w_1=\overline{u}_1-\underline{u}_1$ and
$w_2=\overline{v}_1-\underline{v}_1$. By (A3), (LL), \eqref{Syeq3}, and
the quasimonotone property, we have
\begin{align*}
\frac{\partial w_1}{\partial t}-L_1w_1 +\delta w_1
&=[f_1(x,t,\overline{u},\overline{v}
)+ \delta\overline{u}]-[f_1(x,t,\underline{u},\underline{v}
)+\delta\underline{u}]\\&=[\delta(\overline{u}-\underline{u})+f_1(x,t,\overline{u},\overline{v}
)- f_1(x,t,\overline{u},\underline{v})]\\
&\quad +[f_1(x,t,\overline{u},\underline{v} )-
f_1(x,t,\underline{u},\underline{v} )]\\&\geq 0\;\; \text{in }
\Omega\times {\mathbb{R}},
\end{align*}
and
\begin{align*}
 \mathcal{B}_1 w_1 +\delta w_1
&=[g_1(x,t,\overline{u},\overline{v}
)+\delta\overline{u}]-[g_1(x,t,\underline{u},\underline{v}
)+\delta\underline{u}]\\
&=[\delta (\overline{u}-\underline{u})+g_1(x,t,\overline{u},\overline{v}
)- g_1(x,t,\overline{u},\underline{v})]\\
&+[g_1(x,t,\overline{u},\underline{v} )-
g_1(x,t,\underline{u},\underline{v} )]\\
&\geq 0\quad \text{on }\partial \Omega\times{\mathbb{R}}.
\end{align*}
Since $\sup_{\Omega\times {\mathbb{R}}}\{|w_1(x,t)|\}<\infty$,
it follows from Corollary \ref{Sycor2} that $w_1\geq 0$;
 that is, $\underline{u}_1\leq \overline{u}_1$.
 In a similar way, we prove that $\underline{v}_1\leq\overline{v}_1$.

For $n\geq 2$, a similar argument shows that depending on the choice
of $(u_0,v_0)$, problem \eqref{Syeq3} has solution  either
$(\underline{u}_n, \underline{v}_n)$ or
$(\underline{u}_n, \underline{v}_n) \mathrm{W}^{2,1}_{p,\rm{loc}}(\Omega\times
{\mathbb{R}})\cap L^{\infty}(\Omega \times {\mathbb{R}})\times
\mathrm{W}^{2,1}_{p,\rm{loc}}(\Omega\times
{\mathbb{R}})\cap L^{\infty}(\Omega \times {\mathbb{R}})$
which is such that
$(\underline{u}_{n-1},\underline{v}_{n-1})\leq (u_n,v_n)
\leq(\overline{u}_{n-1}, \overline{v}_{n-1})$.

Indeed, assume by induction that for some $n\geq 2$,
\begin{gather*}
\underline{u}_{n-1}\leq \underline{u}_{n}\leq \overline{u}_{n}\leq
\overline{u}_{n-1},\\
\underline{v}_{n-1}\leq \underline{v}_{n}\leq \overline{v}_{n}\leq
\overline{v}_{n-1}
\end{gather*}
Then, by (A3), (LL) and the quasimonotonicity property, the
functions $w_1=\overline{u}_{n+1}-\underline{u}_{n+1}$ and
$w_2=\overline{v}_{n+1}-\underline{v}_{n+1}$ satisfy
\begin{gather*}
 \frac{\partial w_1}{\partial t}-L_1w_1
+ \delta w_1=[f_1(x,t,\overline{u}_{n},\overline{v}_{n}
)+\delta\overline{u}_{n}]-[f_1(x,t,\underline{u}_n,\underline{v}_n
)+\delta\underline{u}_n]\geq 0\\ \text{in } \Omega\times
{\mathbb{R}},
\\
 \mathcal{B}_1 w_1 +\delta w_1
 =[g_1(x,t,\overline{u}_{n},\overline{v}_{n}
)+\delta \overline{u}_{n}]-[g_1(x,t,\underline{u}_n,\underline{v}_n
)+\delta\underline{u}_n]\geq 0\quad \text{on }
\partial \Omega\times {\mathbb{R}},
\\
\sup_{\Omega\times {\mathbb{R}}}\{|w_1(x,t)|\}<\infty;
\\
 \frac{\partial w_2}{\partial t}-L_2w_2
+\delta w_2=[f_2(x,t,\overline{u}_{n},\overline{v}_{n}
)+\delta\overline{v}_{n}]-[f_2(x,t,\underline{u}_n,\underline{v}_n
)+\delta\underline{v}_n]\geq 0\\
 \text{in } \Omega\times{\mathbb{R}},
\\
 \mathcal{B}_2 w_2 +\delta w_2
 =[g_2(x,t,\overline{u}_{n},\overline{v}_{n}
)+\delta\overline{v}_{n}]-[g_2(x,t,\underline{u}_n,\underline{v}_n
)+\delta\underline{v}_n]\geq 0\quad \text{on }\partial \Omega\times
{\mathbb{R}},
\\
 \sup_{\Omega\times
{\mathbb{R}}}\{|w_2(x,t)|\}<\infty.
\end{gather*}
Using Corollary \ref{Sycor2} we get that $w_i\geq 0$ $(i=1,2)$; that is,
$\underline{u}_{n+1}\leq \overline{u}_{n+1}$  and
$\underline{v}_{n+1}\leq \overline{v}_{n+1}$. Using a similar
argument as above we have that
$\underline{u}_{n}\leq \underline{u}_{n+1}$,
$\overline{u}_{n+1}\leq \overline{u}_{n}$,
$\underline{v}_{n}\leq \underline{v}_{n+1}$, and
$\overline{v}_{n+1}\leq \overline{v}_{n}$.
Thus,
\begin{gather*}
\underline{u}=\underline u_0\leq \underline{u}_1\leq \dots
\leq \underline{u}_n\leq \dots\leq \overline{u}_n\leq \dots
\leq \overline{u}_1\leq \overline{u}_0=\overline{u}
\\
\underline{v}=\underline v_0\leq \underline{v}_1\leq \dots
\leq \underline{v}_n\leq \dots\leq \overline{v}_n\leq \dots\leq
\overline{v}_1\leq \overline{v}_0=\overline{v}
\end{gather*}

Since the sequences $\{u_n\}$ and $\{v_n\}$ (where are $u_n$ represents
either $\underline{u}_n$ or $\overline{u}_n$ and $v_n$ represents either
$\underline{v}_n$ or $\overline{v}_n$) are
monotone then  the pointwise limits
$$
{u^*}(x,t)=\lim_{n\to\infty} u_n(x,t)\quad \text{and}\quad
{v^*}(x,t)= \lim_{n\to\infty} v_n(x,t)
$$
both exist and $\underline{u}\leq {u^*} \leq \overline{u}$,
and $\underline{v}\leq {v^*}\leq \overline{v}$.
We now proceed to show that $ ({u^*},{v^*})$ is  a solution of
 \eqref{Syeq2}.
For that purpose, consider $Q_1=\Omega\times(-1,1)$ and
$Q_2= \Omega\times (-2,2)$. For each
$n\in {\mathbb{N}}$, define $z_n(x,t)= \zeta(t) u_n(x,t)$,
$w_n(x,t)= \zeta(t) v_n(x,t)$, for all
$(x, t)\in \overline{\Omega}\times[-2, 2]$, where
$\zeta\in C^{\infty}({\mathbb{R}}), 0\leq \zeta\leq 1$ and
$\zeta(s)=0$ if $s\leq -2$, $\zeta(s)=1$ if $s\geq -(2-\delta)$ with
$0<\delta<1$. Observe that
$z_n=u_n $  and $w_n= v_n$, in $\overline{\Omega}\times [-1, 1]$,
and satisfy the  linear uncoupled system
\begin{equation}\label{Syeq4}
\begin{gathered}
 \frac{\partial {z}_n}{\partial t}-L_1{z}_n +\delta{z}_n= \frac
{d\zeta}{dt}u_n+\zeta F_{1n}\quad\text{in } \Omega\times
(-2, 2],\\
\frac{\partial {w}_n}{\partial t}-L_2{w}_n +\delta{w}_n=\frac
{d\zeta}{dt}v_n +\zeta F_{2n}\quad\text{in } \Omega\times
(-2, 2],\\
\mathcal{B}_1 {z}_n+\delta{z}_n = \zeta G_{1n}\quad\text{on }
\partial \Omega\times
(-2, 2],\\
\mathcal{B}_2 {w}_n +\delta{w}_n=\zeta G_{2n}\quad\text{on }
\partial \Omega\times
(-2, 2],\\
z_n(x,-2)=0 \quad\text{in } \overline{\Omega},\\
 w_n(x,-2)=0 \quad\text{in } \overline{\Omega},\\
 \sup_{\Omega\times
{\mathbb{R}}}\{|{z}_n(x,t)|, |{w}_n(x,t)|\}<\infty,
\end{gathered}
\end{equation}
where
\begin{gather*}
F_{1n}=f_1(x,t,{u}_{n-1},{v}_{n-1})+\delta{u}_{n-1},\quad
F_{2n}=f_2(x,t,{u}_{n-1},{v}_{n-1})+\delta{v}_{n-1},
\\
G_{1n}=g_1(x,t,{u}_{n-1},{v}_{n-1})+\delta{u}_{n-1},\quad 
G_{2n}=g_2(x,t,{u}_{n-1},{v}_{n-1})+\delta{v}_{n-1}.
\end{gather*}
By the solvability results on linear IBVPs with smooth coefficients
\cite[pp. 341-343]{LSU68S}, it follows that the linear problem \eqref{Syeq4}
have a unique solution
$(z_n, w_n)\in \mathrm{W}_p^{2,  1}(Q_2)\times \mathrm{W}_p^{2,  1}(Q_2) $
(with $p=\frac{N+2}{1-\mu}$). Moreover
\begin{gather}\label{SylPsi1}
|z_n|_{W_p^{2,  1}(Q_2)}
\leq K_0\Big(|\frac{d\zeta}{dt}u_n+\zeta F_{1n}|_{L^p(Q_2)}+|\zeta
G_{1n}|_{W_p^{2-\epsilon-\frac{1}{p}, (2-\epsilon-\frac{1}{p})/2}({\partial
\Omega\times(-2, 2)})}\Big)
\\
\label{SylPsi2}
|w_n|_{W_p^{2,  1}(Q_2)}\leq K_0\Big(|\frac
{d\zeta}{dt}v_n+\zeta F_{2n}|_{L^p(Q_2)}+|\zeta
G_{2n}|_{W_p^{2-\epsilon-\frac{1}{p}, (2-\epsilon-\frac{1}{p})/2}({\partial
\Omega\times(-2, 2)})}\Big),
 \end{gather}
for all $ n\in {\mathbb{N}}$,
 where $K_0$ is a constant which depends on $Q_2$.
Set $V_n=(z_n,w_n)$  with
\[
|V_n|_{W_p^{2,1}(Q_2)}=|z_n|_{W_p^{2,  1}(Q_2)}+|w_n|_{W_p^{2,  1}(Q_2)}.
\]
Observe that for $\epsilon=0$, we get immediately that
$|V_n|_{W_p^{2,  1}(Q_2)}\leq { C}$, for all $n$, since $\varphi_0$
does not depend on $n$. To show that $|V_n|_{W_p^{2,  1}(Q_2)}\leq {C}$
for all $n$ for $\epsilon=1$, we proceed as follows. Using
assumptions (A2) we compute
 $|\zeta G_{i\,n}|_{W_p^{1-\frac{1}{p}, (1-\frac{1}{p})/2}({\partial
\Omega\times(-2, 2)})}$ ($i=1,2$) to get that
\begin{equation}\label{SylPsi3A}
|\zeta G_{i\,n}|_{W_p^{1-\frac{1}{p}, (1-\frac{1}{p})/2}({\partial
\Omega\times(-2, 2)})}\leq
\hat{C}\Big(1+|V_{n-1}|_{W_p^{1-\frac{1}{p}, (1-\frac{1}{p})/2}({\partial
\Omega\times(-2, 2)})}\Big),
\end{equation}
where $\hat{C}$ is independent of $n$ since
$|\zeta G_{i\,n}|_{L^p({\partial \Omega\times(-2, 2)})}
\leq \text{const}$  for all $n \in {\mathbb{N}}$.
Combining \eqref{SylPsi1}, \eqref{SylPsi2}, \eqref{SylPsi3A} we
obtain
$$
|V_n|_{W_p^{2,  1}(Q_2)}\leq
\tilde{C}\Big(1+|V_{n-1}|_{W_p^{1-\frac{1}{p}, (1-\frac{1}{p})/2}({\partial
\Omega\times(-2, 2)})}\Big),
$$
where $\tilde{C}$ is  independent of $n$, but depends on
$|\frac {d\zeta}{dt}v_n+\zeta F_{in}|_{L^p(Q_2)}$,
$|\zeta G_{in}|_{L^p(\partial\Omega \times (-2,2))}$ and the set
$\overline \Omega\times[-2, 2]$.  Using the continuity of the trace operator,
we deduce that
 \begin{equation}\label{SylPsi4}
|V_n|_{W_p^{2,  1}(Q_2)}\leq K\Big(1+|V_{n-1}|_{W_p^{1, 1/2}({
\Omega\times(-2, 2)})}\Big),
\end{equation}
where $K $ does not depend on $n$.
By  the interpolation inequality \eqref{Syintp2}, we get that
\begin{equation}\label{SylPsi5}
|V_n|_{W_p^{2,  1}(Q_2)}
\leq K \Big(1+C\varepsilon |V_{n-1}|_{W_p^{2, 1}(
Q_2)}+\frac{C}{4\varepsilon}|V_{n-1}|_{L^p(Q_2)}\Big)
\end{equation}
From \eqref{SylPsi4} we deduce that
\begin{equation}\label{SylPsi6}
|V_1|_{W_p^{2,  1}(Q_2)}\leq K \Big(1+|\zeta
\overline{V}|_{W_p^{1, 1/2}({ \Omega\times(-2, 2)})}\Big),
\end{equation}
where $\overline{V}$ is either $(\underline{u},\underline{v})$ or
$(\overline{u},\overline{v})$. Combining \eqref{SylPsi5} with
\eqref{SylPsi6} we get
\begin{align*}
|V_2|_{W_p^{2,  1}(Q_2)}
& \leq K\Big(1+C\varepsilon |V_{1}|_{W_p^{2, 1}(
X)}+\frac{C}{4\varepsilon}|V_{1}|_{L^p(Q_2)}\Big)\\
&\leq K\Big(1+KC\varepsilon +KC\varepsilon|\zeta
\overline{V}|_{W_p^{1, 1/2}(Q_2)}+\frac{C}{4\varepsilon}|V_{1}|_{L^p(Q_2)}\Big)
\end{align*}
Proceeding by induction we have  that for every $n\in {\mathbb{N}}$
with $n\geq 2$,
\begin{align*}
|V_n|_{W_p^{2,  1}(Q_2)}
&\leq K \Big(\sum_{i=0}^{n-1}
{(KC\varepsilon)^i}+{(KC\varepsilon)^{n-1}}|\zeta
\overline{V}|_{W_p^{1-\frac{1}{p}, (1-\frac{1}{p})/2}(\partial
\Omega\times(-2, 2))}\\
&\quad + \frac{MC}{4\varepsilon}\sum_{i=0}^
{n-2}{(KC\varepsilon)^i}\Big),
\end{align*}
where  $K$  is independent of $n$,  and the constant
$M\geq|V_{n}|_{L^p(Q_2)}$ for all $n\in {\mathbb{N}}$.
Therefore, we obtain the following estimate which involves a geometric series
$$
|V_n|_{W_p^{2,  1}(Q_2)}\leq \Big(K+ K|\zeta
\overline{V}|_{W_p^{1-\frac{1}{p}, (1-\frac{1}{p})/2}(\partial
\Omega\times(-2, 2))}+ \frac{MCK}{4\varepsilon}\Big)\sum_{i=0}^
{\infty}(KC\varepsilon)^i.
$$
Thus,
$|V_n|_{W_p^{2,  1}(Q_2)}\leq {C} $ for all $n\in {\mathbb{N}}$,
provided $\varepsilon$ is chosen sufficiently small such that
$KC\varepsilon<1$.
It follows that
$|z_n|_{W_p^{2,  1}(Q_2)}\leq {C}, \,|w_n|_{W_p^{2, 1}(Q_2)}\leq {C}$.

 Now, we need to show that in $Q_1$, the sequence
$\{V_n\}=\{(u_n,v_n)\}$ has a
subsequence which converges to a solution of  problem \eqref{Syeq2}.
Indeed, define $T: \big(W_p^{2,  1}(Q_1),|.|_{W_p^{2,
1}(Q_1)}\big) \to (L^p(Q_1),|.|_{L^p(Q_1)})$ by
$$
T(v)=\frac{\partial v}{\partial t}-Lv+\delta v.
$$
Hence, $T$ is (weakly) closed.  Since $W_p^{2,  1}(Q_1)$ is a
reflexive space  and $|V_n|_{W_p^{2,  1}(Q_1)}\leq \tilde{C} $ for
all $n$, there exist subsequences $\{u_n\}$ and  $\{v_n\}$ such that
$$
u_n\rightharpoonup \tilde{u}_1\quad \text{and} \quad
v_n\rightharpoonup \tilde{v}_1\quad \text{in } W_p^{2,  1}(Q_1).
$$
By the compact embedding of $ W_p^{2,  1}(Q_1)$ into
$C^{1+\mu, (1+\mu)/2}(\overline{Q}_1)$, it follows that there exist
 subsequences $\{u_{1n}\}$ and $\{v_{1n'}\}$ such that
$u_{1n}\to \tilde{u}_1$, $v_{1n'}\to \tilde{v}_1$ in
$C^{1+\mu, (1+\mu)/2}(\overline{Q}_1)$. Since $v_{n}$ is a
monotone sequence, it follows that the  sequence $v_{n}\to
\tilde{v}_1$ uniformly in $\overline{Q}_1$. Therefore, any
subsequence of $\{v_{n}\}$  converges uniformly in $\overline{Q}_1$
to $\tilde{v}_1$  in particular  $\{v_{1n}\}$.
Moreover, since $T$ is (weakly) closed and
$T(u_{1n})=F_{1n}\to f_1(.,.,\tilde{u}_1 ,\tilde{v}_1)+k\,
\tilde{u}_{1}$
uniformly in $\overline{Q}_1$, it follows that
$T(\tilde{u}_1)=f_1(.,., \tilde{u}_1 ,\tilde{v}_1)+k\, \tilde{u}_{1}$.
In addition,
$\mathcal{B}_1 {u_{1n}}+\epsilon k\,u_{1n}\to \mathcal{B}_{1}
\tilde{u}_1+\epsilon k\, \tilde{u}_{1}$ in
$C^{\mu, \mu/2}(\bar{\partial \Omega\times[-1, 1]})$, and
$\mathcal{B}_1 {u_{1n}}+\epsilon k \,u_{1n}=G_{1n} \to
g_1(.,., \tilde{u}_1 ,\tilde{v}_1)+k\, \tilde{u}_{1}$ uniformly on
$\partial \Omega\times[-1, 1]$, we get
$\mathcal{B}_1\tilde{u}_1+\epsilon k\,\tilde{u}_{1}=g_1(.,., \tilde{u}_1
,\tilde{v}_1)+k\, \tilde{u}_{1}$. Therefore,
$\tilde{u}_1$ satisfies the following problem
\begin{gather*}
\frac{\partial \tilde{u}_1}{\partial t}-L\tilde{u}_1+k\tilde{u}_1
=f_1(x,t, \tilde{u}_1 ,\tilde{v}_1)+k\, \tilde{u}_{1}\quad
\text {in }\Omega\times(-1, 1),
\\
 \mathcal{B}_1\tilde{ u}_1+\epsilon
k\,\tilde{u}_{1}=g_1(x,t, \tilde{u}_1 ,\tilde{v}_1)+k\,
\tilde{u}_{1} \quad \text {on }\partial \Omega\times[-1, 1],\\
\sup_{\Omega\times [-1, 1]}|\tilde{u}_1(x,t)|<\infty.
\end{gather*}
Using  similar arguments as above for $v_n$, we obtain
\begin{gather*}
\frac{\partial \tilde{v}_1}{\partial t}-L\tilde{v}_1+k\tilde{v}_1
=f_2(x,t, \tilde{u}_1 ,\tilde{v}_1)+k\, \tilde{v}_{1} \quad
\text {in }\Omega\times(-1,1),\\
\mathcal{B}_1\tilde{ v}_1+\epsilon k\,\tilde{v}_{1}=g_2(x,t,
\tilde{u}_1 ,\tilde{v}_1)+k\, \tilde{v}_{1} \quad
\text{on }\partial \Omega\times[-1, 1],\\
\sup_{\Omega\times [-1, 1]}|\tilde{v}_1(x,t)|<\infty.
\end{gather*}
For $n\geq 2$, let $Q_n=\Omega\times(-n, n)$.
Consider the subsequence $\{(u_{(n-1)k}, v_{(n-1)k})\}$
and use similar arguments to the above to extract a subsequence
$\{(u_{n\, k}, v_{n\,k})\}$ of $\{(u_{(n-1)k}, v_{(n-1)k})\}$ such that
it converges to  $(\tilde{u}_n,\tilde{v}_n)$ in
$C^{1+\mu, (1+\mu)/2}(\overline{\Omega}\times [-n,n])
\times C^{1+\mu, (1+\mu)/2}(\overline{\Omega}\times [-n,n])$ which
satisfies
\begin{gather*}
\frac{\partial \tilde{u}_n}{\partial t}-L\tilde{u}_n+k\tilde{u}_n
=f_1(x,t, \tilde{u}_n ,\tilde{v}_n)+k \tilde{u}_{n} \quad\text {in }
 \Omega\times(-n, n),\\
\frac{\partial \tilde{v}_n}{\partial t}-L\tilde{v}_n+k\tilde{v}_n
=f_2(x,t,\tilde{u}_n ,\tilde{v}_n)+k\, \tilde{v}_{n} \quad
\text {in }\Omega\times(-n, n),\\
\mathcal{B}_1\tilde{ u}_n+\epsilon k\,\tilde{u}_{n}
=g_2(x,t,\tilde{u}_n ,\tilde{v}_n)+k \tilde{u}_{n} \quad\text{on }
\partial \Omega\times[-n, n],\\
 \mathcal{B}_1\tilde{ v}_n+\epsilon k\,\tilde{v}_{n}
=g_2(x,t, \tilde{u}_n ,\tilde{v}_n)+k\tilde{v}_{n} \quad\text{on }
\partial \Omega\times[-n, n],\\
\sup_{\Omega\times [-n, n]}|\tilde{u}_n(x,t),\tilde{v}_n(x,t)|<\infty.
\end{gather*}
Note that by  construction,
$(u_n, v_n)|_{\Omega\times [-(n-1), n-1]}=(u_{n-1},v_{n-1})$ for all
$n\ge 2$; that is, $(u_n, v_n)$ is an extension of $(u_{n-1},v_{n-1})$.
Now, set $ V_n=\{(u_n,v_n)\}$. By  the diagonalization argument,
choose the sequence $\{V_{jj}\}$ located on the `diagonal'.
Observe that  $V_{jj}\in \{V_{nk}\}=\{(u_{nk},v_{nk})\}$ for every $n\leq j$,
and hence $\{V_{jj}\}$ is a subsequence of $\{V_n\}$. We shall prove that
the sequence $\{V_{jj}\}$ converges to a solution
$V^*$ of problem \eqref{Syeq2}.
Indeed, let   $\overline{\Omega}\times [-n,n] $ and $\varepsilon>0 $.
Since $\{V_{nk}\}$ converges to $\tilde{V}_n= (\tilde{u}_n,\tilde{v}_n)$ in
$C^{1+\mu, (1+\mu)/2}(\overline{\Omega}\times [-n,n])$, there exists
$ N\in{\mathbb{N}}$ such that for all
$k\geq N,  |V_{nk}-\tilde{V}_k|_{C^{1+\mu, (1+\mu)/2}(\overline{\Omega}\times
[-n,n])}  <\varepsilon$. Using the fact that $V_{jj}\in \{V_{nk}\}$
 for all $j\geq n$, we get that  for all $j\geq \max\{n, N\}$,
$|V_{jj}-\tilde{V}_n|_{C^{1+\mu, (1+\mu)/2}(\overline{\Omega}\times [-n,n])}
<\varepsilon$.

Thus, $\{V_{jj}\}$ is subsequence of $\{V_n\}$ which  converges
(on every compact set) to a function $\tilde{V}$ in
$C^{1+\mu, (1+\mu)/2}(\overline{\Omega}\times [-n,n])$,
where $\tilde{V}|_{\overline{\Omega}\times [-n,n]}=\tilde{V}_n$, so that
$\tilde{V}\in C_{\rm loc}^{1+\mu, (1+\mu)/2}(\overline{\Omega}\times
{\mathbb{R}})\cap W_{p,\text{loc}}^{2, 1}(\Omega\times{\mathbb{R}} )$ and
$\sup_{\Omega\times {\mathbb{R}}}|\tilde{V}|\leq M$.
 Moreover, $\tilde{V}=(\tilde{u},\tilde{v})$ satisfies the problem
\eqref{Syeq2}.  By  uniqueness of the limit we get
$\tilde{V}=V^*$=$({u^*}, {v^*})$.
Thus, $( u^* ,v^*)$ is a solution of problem \eqref{Syeq1} and
$(\underline{u},\underline{v})\leq ( u^* ,v^*)\leq (\overline{u},
\overline{v})$. The proof is complete.
\end{proof}

To prove Theorem \ref{Syt1}, we will use (an improved version of)
a nonlinear approximation argument inspired by the one considered
in \cite{NM11} (see also \cite{HA71}). However, the main difficulty
lies in the obtainment of required a priori estimates since there
is a lack of compactness herein. We will therefore need the preliminary
lemmas that as proved below. For sake of discussion, we will present
the rest of the proof for the quasimonotone nonincreasing case.

\begin{lemma}\label{transfo1b}
Let $f_1,f_2$ satisfy {\rm{(A1)}}; that is,
$f_1,\,f_2\in L^\infty_{\rm{car}}(\Omega\times {\mathbb{R}}\times\mathbb{R})$.
Then, for every $r\in\mathbb{R}$ with $r>0$, there are  continuous
functions $m_1, m_2:[-r,r]\times[-r,r]\to\mathbb{R}$ such that $m_i(\cdot,v)$
is nondecreasing on $[-r,r]$, $m_i(u,\cdot)$ is nonincreasing on
$[-r,r]$, $m_i(u,v)=-m_i(v,u)$ on $[-r,r]\times[-r,r]$, and
\begin{equation}\label{transfo2}
\sup_{\Omega\times\mathbb{R}}|f_1(x,t,u,w)-f_1(x,t,v,w)|\le m_1(u,v)
\end{equation}
for all $u,v\in[-r,r]$ with $u\ge v$, and
\begin{equation}\label{transfo21}
\sup_{\Omega\times\mathbb{R}}|f_2(x,t,w,u)-f_2(x,t,w,v)|\le m_2(u,v)
\end{equation}
for all $u,v\in[-r,r]$ with $u\ge v$.
\end{lemma}

The proof of Lemma \ref{transfo1} is similar to \cite[Lemma 3.4]{NM11}.
Setting (for instance)
$$
r=\max\big(|\underline{u}|_{_{L^\infty(\Omega\times\mathbb{R})}},
|\overline{u}|_{_{L^\infty(\Omega\times\mathbb{R})}},
\underline{u}|_{_{L^\infty(\Omega\times\mathbb{R})}},
\overline{u}|_{_{L^\infty(\Omega\times\mathbb{R})}}\big)+2,
$$
it follows from the Stone-Weierstrass Approximation Theorem that for
every $n\in\mathbb{N}$ there is a Lipschitz continuous function
$m_{i,n}:[-r,r]\times[-r,r]\to\mathbb{R}$ such that
\begin{equation}\label{approx1}
|m_i(u,v)-m_{i,n}(u,v)|<\frac1n
\end{equation}
for all $(u,v)\in[-r,r]\times[-r,r]$.

Now, consider the modified  problems
\begin{equation}\label{NLmp2}
\begin{gathered}
\frac{\partial u}{\partial t}(x,t)-Lu (x,t)= {f}_1(x,t,\underline{u},
\underline{v})+m_1(\underline{u}, u) \quad \text{a.e. in } \Omega\times
{\mathbb{R}},\\
\frac{\partial v}{\partial t}(x,t)-Lv (x,t)= {f}_2(x,t,\underline{u},
\underline{v})+m_2(\underline{v}, v) \quad \text{a.e. in } \Omega\times
{\mathbb{R}},\\
\mathcal{B}_\epsilon u = g_1(x,t,\underline{u},
\underline{v})+\rho_1(\underline{u}- u) \quad \text{on }
\partial \Omega\times {\mathbb{R}},\\
\mathcal{B}_\epsilon v = g_2(x,t,\underline{u}, \underline{v})
+\rho_2(\underline{v}- v) \quad \text{on } \partial \Omega\times {\mathbb{R}},\\
 \sup_{\Omega\times{\mathbb{R}}}\{|u(x,t)|,|v(x,t)| \}<\infty,
\end{gathered}
\end{equation}
and
\begin{equation}\label{NLmp1}
\begin{gathered}
 \frac{\partial u}{\partial t}(x,t)-Lu (x,t)= {f}_1(x,t,\overline{u},
\overline{v})+m_1(\overline{u}, u) \quad \text{a.e. in } \Omega\times
{\mathbb{R}},\\
\frac{\partial v}{\partial t}(x,t)-Lv (x,t)= {f}_2(x,t,\overline{u},
\overline{v})+m_2(\overline{v}, v) \quad \text{a.e. in } \Omega\times
{\mathbb{R}},\\
\mathcal{B}_1 u = g_1(x,t,\overline{u}, \overline{v})+\varrho_1(\overline{u}- u)
\quad \text{on } \partial \Omega\times {\mathbb{R}},\\
\mathcal{B}_2 v = g_2(x,t,\overline{u}, \overline{v})+\varrho_2(\overline{v}- v)
\quad \text{on } \partial \Omega\times {\mathbb{R}},\\
 \sup_{\Omega\times {\mathbb{R}}}\{|u(x,t)|,|v(x,t)| \}<\infty.
\end{gathered}
\end{equation}
Define the functions $\hat{f}_i, \check{f}_i\in L^\infty_{\rm{car}}(\Omega\times
{\mathbb{R}}\times [\underline{u},\overline{u}])$ $(i=1,2)$  by
\begin{gather*}
\hat{f}_1(x,t, u,\underline{v}):= {f}_1(x,t,\underline{u}, \underline{v})
+m_1(\underline{u}, u);
\\
\hat{f}_2(x,t, \underline{u}, {v}):= {f}_2(x,t,\underline{u},
\underline{v})+m_2(\underline{v}, v)
\\
\check{f}_1(x,t, u,\overline{v}):={f}_1(x,t,\overline{u},
\overline{v})+m_1(\overline{u}, u);
\\
\check{f}_2(x,t, \overline{u},v):={f}_2(x,t,\overline{u},
\overline{v})+m_2(\overline{v}, v).
\end{gather*}
Define the functions $\hat{g}_i\;\check{g}_i$, which satisfy  condition
(A2), by
\begin{gather*}
\hat{g}_1(x,t, u,\underline{v}):= g_1(x,t,\underline{u},
\underline{v})+\rho_1(\underline{u}- u);
\\
\hat{g}_2(x,t, \underline{u}, {v}):= g_2(x,t,\underline{u},
\underline{v})+\rho_2(\underline{v}- v)
\\
\check{g}_1(x,t, u,\overline{v}):=g_1(x,t,\overline{u}, \overline{v})
+\varrho_1(\overline{u}- u);
\\
\check{g}_2(x,t, \overline{u},v):=g_2(x,t,\overline{u}, \overline{v})
+\varrho_2(\overline{v}- v).
\end{gather*}
Observe that $\hat{f}_1, \hat{g}_1$ ($\hat{f}_2, \hat{g}_2$)
are nonincreasing in $u\in(-\infty, \underline{u}]$
(in $v\in (-\infty, \underline{v}]$), and
$\check{f}_1,\check{g}_1$ ($\check{f}_2,\check{g}_2$) are nonincreasing
in $u\in [\overline{u},\infty)$  (in $v\in [\underline{v},\infty)$).
Moreover, by using Lemma \ref{transfo1}, \eqref{transfo2} and
\eqref{transfo21}, they satisfy the following inequalities:
\begin{gather*}
\hat{f}_1(x,t, \cdot,\underline{v})\leq {f}_1(x,t, \cdot,\underline{v})
\leq \check{f}_1(x,t, \cdot,\overline{v});
\\
\hat{g}_1(x,t, \cdot,\underline{v})\leq {g}_1(x,t, \cdot,\underline{v})
\leq \check{g}_1(x,t, \cdot,\overline{v})
\end{gather*}
on $[\underline{u},\overline{u}]$.
\begin{gather*}
\hat{f}_2(x,t, \underline{u}, \cdot)\leq \hat{f}_2(x,t, \underline{u}, \cdot)
\leq \check{f}_2(x,t, \overline{u},\cdot);
\\
\hat{g}_2(x,t, \underline{u}, \cdot)\leq \hat{g}_2(x,t, \underline{u}, \cdot)
\leq \check{g}_2(x,t, \overline{u},\cdot)
\end{gather*}
on $[\underline{v},\overline{v}]$.

We will show that problem \eqref{NLmp2} and problem \eqref{NLmp1} have unique
solutions in $[(\underline{u},\underline{v}),(\overline{u},\overline{v})]$.
In order to accomplish this, we first need the following lemma.

\begin{lemma}\label{lemdelta}
Assume that {\rm (A1)--(A3)} are satisfied and that
$(\underline{u},\underline{v})$
and $(\overline{u},\overline{v})$ are subsolution and supersolution of
problem \eqref{Syeq1} with $\underline{u}\leq \overline{u}$ and
$\underline{v}\leq \overline{v}$ . Let
$\delta>0$ and $\overline{u}_\delta:=\overline{u}+\delta z$,
$\underline{u}_\delta:=\underline{u}-\delta z$,
$\overline{v}_\delta:=\overline{v}+\delta z$,
$\underline{v}_\delta:=\underline{v}-\delta z$, where $z$ is the
(unique) solution of the linear boundary value problem
\begin{equation}\label{z}
\frac{\partial u}{\partial t}(x,t)-Lu (x,t)=1 \quad\text{a.e in}\quad \Omega\times
{\mathbb{R}}, \; \mathcal{B} u =1+\epsilon
\quad\text{on }\partial \Omega\times {\mathbb{R}},
\end{equation}
 and $\sup_{\Omega\times {\mathbb{R}}}|u(x,t)|<\infty$. Then the boundary
value problems
\begin{equation}\label{NLmp4}
\begin{gathered}
\frac{\partial u}{\partial t}(x,t)-Lu (x,t)=
 f(x,t,\underline{u},\underline{v})+m_1(\underline{u}_\delta,u)
\quad \text{a.e. in }\Omega\times{\mathbb{R}},\\
 \frac{\partial v}{\partial t}(x,t)-Lv (x,t)=
 f(x,t,\underline{u},\underline{v})+m_2(\underline{v}_\delta,v)  \quad
\text{a.e. in } \Omega\times {\mathbb{R}},\\
\mathcal{B}_1 u = g_1(x,t,\underline{u}, \underline{v})
+\varrho_1(\underline{u}_\delta-u) \quad
\text{on } \partial \Omega\times {\mathbb{R}},\\
\mathcal{B}_2 v = g_2(x,t,\underline{u}, \underline{v})
+\varrho_2(\underline{v}_\delta-v) \quad
\text{on } \partial \Omega\times {\mathbb{R}},\\
\sup_{\Omega\times {\mathbb{R}}}\{|u(x,t)|,|v(x,t)|\}<\infty,
\end{gathered}
\end{equation}
and
\begin{equation}\label{NLmp3}
\begin{gathered}
\frac{\partial u}{\partial t}(x,t)-Lu (x,t)=
 f(x,t,\overline{u},\overline{v})+m_1(\overline{u}_\delta,u)  \quad
 \text{a.e. in } \Omega\times{\mathbb{R}},\\
 \frac{\partial v}{\partial t}(x,t)-Lv (x,t)=
 f(x,t,\overline{u},\overline{v})+m_2(\overline{v}_\delta,v)  \quad
 \text{a.e. in } \Omega\times{\mathbb{R}},\\
\mathcal{B}_1 u = g_1(x,t,\overline{u}, \overline{v})
+\varrho_1(\overline{u}_\delta-u) \quad
\text{on } \partial \Omega\times {\mathbb{R}},\\
\mathcal{B}_2 v = g_2(x,t,\overline{u}, \overline{v})
+\varrho_2(\overline{v}_\delta-v) \quad
 \text{on } \partial \Omega\times {\mathbb{R}},\\
 \sup_{\Omega\times
{\mathbb{R}}}\{|u(x,t)|,|v(x,t)|\}<\infty,
\end{gathered}
\end{equation}
have unique solutions $(\hat{u}_\delta,\hat{v}_\delta)$ and
 $(\check{u}_\delta,\check{v}_\delta)$
 respectively such that
 $$
\underline{u}_\delta\leq\hat{u}_\delta\le
\check{u}_\delta \leq\overline{u}_\delta \quad\text{and}\quad
\underline{v}_\delta\leq\hat{v}_\delta\le
\check{v}_\delta \leq\overline{v}_\delta.
$$
\end{lemma}

\begin{proof}
Let us define the functions
\begin{gather*}
\hat{f}_{1\delta}(x,t, u,\underline{v}):= {f}_1(x,t,\underline{u},
\underline{v})+m_1(\underline{u}_{\delta}, u);
\\
\hat{f}_{2\delta}(x,t, \underline{u}, {v}):= {f}_2(x,t,\underline{u},
\underline{v})+m_2(\underline{v}_{\delta}, v)
\\
\check{f}_{1\delta}(x,t, u,\overline{v}):={f}_1(x,t,\overline{u},
\overline{v})+m_1(\overline{u}_{\delta}, u);
\\
\check{f}_{2\delta}(x,t, \overline{u},v):={f}_2(x,t,\overline{u},
\overline{v})+m_2(\overline{v}_{\delta}, v).
\\
\hat{g}_{1\delta}(x,t, u,\underline{v}):= g_1(x,t,\underline{u}, \underline{v})
+\rho_1(\underline{u}_{\delta}- u);
\\
\hat{g}_{2\delta}(x,t, \underline{u}, {v}):= g_2(x,t,\underline{u},
\underline{v})+\rho_2(\underline{v}_{\delta}- v);
\\
\check{g}_{1\delta}(x,t, u,\overline{v}):=g_1(x,t,\overline{u},
\overline{v})+\varrho_1(\overline{u}_{\delta}- u) ;
\\
\check{g}_{2\delta}(x,t, \overline{u},v):=g_2(x,t,\overline{u}, \overline{v})
+\varrho_2(\overline{v}_{\delta}- v).
\end{gather*}

From the monotonicity properties of $m_i$, one has that
$\hat{f}_{1\delta},\hat{g}_{1\delta}$ ($\hat{f}_{2\delta},\hat{g}_{2\delta}$)
are nonincreasing functions of $u$ for $u\geq   \underline{u}_{\delta}$
(of $v$ for $v\geq   \underline{v}_{\delta}$), and
$\check{f}_{1\delta}, \check{g}_{1\delta}$
($\check{f}_{2\delta}, \check{g}_{2\delta}$) are nonincreasing functions
of $u$ for $u\leq   \overline{u}_{\delta}$ (of $v$ for
$v\leq   \overline{v}_{\delta}$).

Using the definitions and the monotonicity properties of $m_i$,
the quasimonotonicity property of $f_i$, and the fact that
$\hat{f}_1(x,t, \cdot,\underline{v})\leq \check{f}_1(x,t, \cdot,\overline{v})$
on $[\underline{u}, \overline{u}]$ and
$\hat{f}_2(x,t, \underline{u}, \cdot)\leq \check{f}_2(x,t, \overline{u},\cdot)$
on $[\underline{v}, \overline{v}]$, it follows that
\begin{equation}\label{ineqf1}
\hat{f}_{1\delta}(x,t, u,\underline{v})
\le \check{f}_{1\delta}(x,t, u,\overline{v})\quad\text{for all }
 u\in [\underline{u}_{\delta}, \overline{u}_{\delta}]
 \end{equation}
and
 \begin{equation}\label{ineqf2}
 \hat{f}_{2\delta}(x,t, \underline{u},{v})\le \check{f}_{2\delta}(x,t,
\overline{u},{v})\quad \text{for all }
 v\in [\underline{v}_{\delta}, \overline{v}_{\delta}].
\end{equation}
Now, by using  \eqref{z} and the fact that $(\underline{u},\underline{v})$
and $(\overline{u},\overline{v})$ are sub- and super-solutions,
it is seen that
\begin{equation}\label{strictsubs}
\begin{gathered}
\frac{\partial \underline{u}_{\delta}}{\partial t}(x,t)
-L\underline{u}_{\delta} (x,t)-\hat{f}_{1\delta}(x,t, \underline{u}_{\delta},
\underline{v})\leq -\delta   \quad\text{a.e. in } \Omega\times {\mathbb{R}},\\
\frac{\partial \underline{v}_{\delta} }{\partial t}(x,t)-L\underline{v}_{\delta}
(x,t)-\hat{f}_{2\delta}(x,t, \underline{u}, \underline{v}_{\delta})
\leq -\delta   \quad\text{a.e. in } \Omega\times {\mathbb{R}},
\\
\mathcal{B}_1 \underline{u}_{\delta}-\hat{g}_{1\delta}
(x,t, \underline{u}_{\delta},\underline{v})\leq -2\delta \quad
\text{on } \partial \Omega\times {\mathbb{R}},
\\
\mathcal{B}_2 \underline{v}_{\delta} -\hat{g}_{2\delta}(x,t, \underline{u},
\underline{v}_{\delta})\leq -2\delta \quad
\text{on } \partial \Omega\times {\mathbb{R}},\\
\sup_{\Omega\times {\mathbb{R}}}\{|\underline{u}_{\delta}(x,t)|,
|\underline{v}_{\delta}(x,t)|\}<\infty,
\end{gathered}
\end{equation}
and
\begin{equation}\label{strictsups}
\begin{gathered}
 \frac{\partial \overline{u}_{\delta}}{\partial t}(x,t)
-L\overline{u}_{\delta} (x,t)-\check{f}_{1\delta}(x,t, \overline{u}_{\delta},
\overline{v})\geq \delta  \quad\text{a.e. in } \Omega\times
{\mathbb{R}},\\
 \frac{\partial \overline{v}_{\delta} }{\partial t}(x,t)-L\overline{v}_{\delta}
 (x,t)-\check{f}_{2\delta}(x,t, \overline{u}, \overline{v}_{\delta})\geq \delta
  \quad\text{a.e. in } \Omega\times {\mathbb{R}},\\
\mathcal{B}_1 \overline{u}_{\delta}-\check{g}_{1\delta}(x,t,
\overline{u}_{\delta},\overline{v})\geq 2\delta \quad
\text{on } \partial \Omega\times {\mathbb{R}},\\
\mathcal{B}_2 \overline{v}_{\delta} -\check{g}_{2\delta}(x,t, \overline{u},
\overline{v}_{\delta})\geq 2\delta
\quad \text{on } \partial \Omega\times {\mathbb{R}},\\
\sup_{\Omega\times {\mathbb{R}}}\{|\overline{u}_{\delta}(x,t)|,
|\overline{v}_{\delta}(x,t)|\}<\infty.
\end{gathered}
\end{equation}
Therefore $(\underline{u}_{\delta},\underline{v}_{\delta})$ is a strict
 subsolution of problem \eqref{NLmp4} and
$(\overline{u}_{\delta},\overline{v}_{\delta})$ is a strict supersolution
of problem \eqref{NLmp3}.
The functions $\hat{f}_1,\hat{f}_2,\check{f}_1, \check{f}_2,\hat{g}_1,
\hat{g}_2,\check{g}_1$, and $\check{g}_2 $ defined on
$ \overline{\Omega}\times {\mathbb{R}}\times
[\underline{u}_{\delta},\overline{u}_{\delta}]
\times[\underline{v}_{\delta}, \overline{v}_{\delta}]$
may be approximated, for $n\in {\mathbb{N}}$, by the following functions
\begin{gather*}
\hat{f}_{1n}(x,t, u,\underline{v}):= {f}_1(x,t,\underline{u}, \underline{v})
+m_{1n}(\underline{u}_{\delta}, u); \\
\hat{f}_{2n}(x,t, \underline{u}, {v}):= {f}_2(x,t,\underline{u},
\underline{v})+m_{2n}(\underline{v}_{\delta}, v)
\\
\check{f}_{1n}(x,t, u,\overline{v}):={f}_1(x,t,\overline{u},
\overline{v})+m_{1n}(\overline{u}_{\delta}, u);\\
\check{f}_{2n}(x,t, \overline{u},v):={f}_2(x,t,\overline{u},
\overline{v})+m_{2n}(\overline{v}_{\delta}, v).
\\
\hat{g}_{1n}(x,t, u,\underline{v}):= g_1(x,t,\underline{u}, \underline{v})
+\rho_{1n}(\underline{u}_{\delta}- u); \\
\hat{g}_{2n}(x,t, \underline{u}, {v}):= g_2(x,t,\underline{u},
\underline{v})+\rho_{2n}(\underline{v}_{\delta}- v)
\\
\check{g}_{1n}(x,t, u,\overline{v}):=g_1(x,t,\overline{u}, \overline{v})
+\varrho_{1n}(\overline{u}_{\delta}- u);\\
\check{g}_{2n}(x,t, \overline{u},v):=g_2(x,t,\overline{u}, \overline{v})
+\varrho_{2n}(\overline{v}_{\delta}- v),
\end{gather*}
where $m_{in}$ are the Lipschitz approximation of $m_i$ satisfying
\eqref{approx1}. Moreover, it is easy to check that
 $(\underline{u}_{\delta}, \underline{v}_{\delta})$ and
$(\overline{u}_{\delta},\overline{v}_{\delta})$ are still supersolution
and subsolution of the following approximating equations
\begin{equation}\label{NLmp41}
\begin{gathered}
\frac{\partial u}{\partial t}(x,t)-Lu (x,t)=
 f_1(x,t,\underline{u},\underline{v})+m_{1n}(\underline{u}_\delta,u)
 \quad\text{a.e. in } \Omega\times {\mathbb{R}},\\
 \frac{\partial v}{\partial t}(x,t)-Lv (x,t)=
 f_2(x,t,\underline{u},\underline{v})+m_{2n}(\underline{v}_\delta,v)
\quad\text{a.e. in } \Omega\times {\mathbb{R}},\\
\mathcal{B}_1 u = g_1(x,t,\underline{u}, \underline{v})
+\varrho_1(\underline{u}_\delta-u)
\quad \text{on } \partial \Omega\times {\mathbb{R}},\\
\mathcal{B}_2 v = g_2(x,t,\underline{u}, \underline{v})
+\varrho_2(\underline{v}_\delta-v)
\quad \text{on } \partial \Omega\times {\mathbb{R}},\\
\sup_{\Omega\times {\mathbb{R}}}\{|u(x,t)|,|v(x,t)|\}<\infty,
\end{gathered}
\end{equation}
and
\begin{equation}\label{NLmp31}
\begin{gathered}
 \frac{\partial u}{\partial t}(x,t)-Lu (x,t)=
 f_1(x,t,\overline{u},\overline{v})+m_{1n}(\overline{u}_\delta,u)
\quad\text{a.e. in } \Omega\times {\mathbb{R}},\\
 \frac{\partial v}{\partial t}(x,t)-Lv (x,t)=
 f_2(x,t,\overline{u},\overline{v})+m_{2n}(\overline{v}_\delta,v)
 \quad\text{a.e. in } \Omega\times{\mathbb{R}},\\
\mathcal{B}_1 u = g_1(x,t,\overline{u}, \overline{v})
+\varrho_1(\overline{u}_\delta-u)
\quad \text{on } \partial \Omega\times {\mathbb{R}},\\
\mathcal{B}_2 v = g_2(x,t,\overline{u}, \overline{v})
+\varrho_2(\overline{v}_\delta-v)
\quad \text{on } \partial \Omega\times {\mathbb{R}},\\
 \sup_{\Omega\times {\mathbb{R}}}\{|u(x,t)|,|v(x,t)|\}<\infty.
\end{gathered}
\end{equation}
Since $(\underline{u}_{\delta},\underline{v}_{\delta})
\leq (\overline{u}_{\delta},\overline{v}_{\delta})$ and the
functions $\hat{f}_{in}$ and $\check{f}_{in}$ satisfy  the
(LL)-condition in
$[\underline{u}_{\delta},\overline{u}_{\delta}]
\times [\underline{v}_{\delta},\overline{v}_{\delta}]$, it follows from
Proposition \ref{Syp11} that there is a solution
$(\hat{u}_{\delta,n},\hat{v}_{\delta,n})$ of problem \eqref{NLmp41}
and a solution $(\check{u}_{\delta n},\check{v}_{\delta n})$
of problem \eqref{NLmp31} such that
$(\underline{u}_{\delta},\underline{v}_{\delta})
\leq (\hat{u}_{\delta n},\hat{v}_{\delta n})
\leq (\overline{u}_{\delta},\overline{v}_{\delta})$
and $(\underline{u}_{\delta},\underline{v}_{\delta})
\leq (\check{u}_{\delta n},\check{v}_{\delta n})
\leq (\overline{u}_{\delta},\overline{v}_{\delta})$.

Now we proceed to show that (a relabeled subsequence of)
$(\hat{u}_{\delta n},\hat{v}_{\delta n})$ converges
(uniformly on compact sets) to a solution
$(\hat{u}_{\delta},\hat{v}_{\delta})$ of problem \eqref{NLmp4} with
$(\underline{u}_{\delta},\underline{v}_{\delta})
\leq (\hat{u}_{\delta},\hat{v}_{\delta}) \leq
(\overline{u}_{\delta},\overline{v}_{\delta})$.

For that purpose, consider $Q_1=\Omega\times (-1, 1)$ and
$Q_2= \Omega\times (-2, 2). $ For each $n\in {\mathbb{N}}$, define
$\hat{z}_n(x,t)= \zeta(t) \hat{u}_{\delta, n}(x,t)$,
$\hat{w}_n(x,t)= \zeta(t) \hat{v}_{\delta, n}(x,t)$, for all
 $(x, t)\in \overline{\Omega}\times [-2,  2]$, where
$\zeta\in C^{\infty}({\mathbb{R}}),  0\leq \zeta\leq 1$ and
$\zeta(s)=0$  if  $s\leq -2$,  $\zeta(s)=1$
if $s\geq -(2-\delta)$ with $0<\delta<1$. Observe that
$z_n=u_n $  and $w_n= v_n$, in $\overline{\Omega}\times [-1,  1]$,
and satisfy the  linear uncoupled system
\begin{equation}\label{Syeq4b}
\begin{gathered}
 \frac{\partial {\hat{z}}_n}{\partial t}-L_1{\hat{z}}_n = \frac
{d\zeta}{dt}\hat{u}_{\delta n}+\zeta \hat{f}_{1n}\quad\text{in } \Omega\times
(-2, 2],\\
\frac{\partial {\hat{w}}_n}{\partial t}-L_2{\hat{w}}_n =\frac
{d\zeta}{dt} \hat{v}_{\delta n} +\zeta \hat{f}_{2n}\quad\text{in } \Omega\times
(-2, 2],\\
\mathcal{B}_1 {\hat{z}}_n = \zeta \hat{g}_{1n}\quad\text{on }
\partial \Omega\times
(-2, 2],\\
\mathcal{B}_2 {\hat{w}}_n =\zeta \hat{g}_{2n}\quad\text{on }
\partial \Omega\times
(-2, 2],\\
\hat{z}_n(x,-2)=0 \quad\text{in } \overline{\Omega},\\
\hat{ w}_n(x,-2)=0 \quad \text{in } \overline{\Omega},\\
 \sup_{\Omega\times {\mathbb{R}}}\{|{\hat{z}}_n(x,t)|,
|{\hat{w}}_n(x,t)|\}<\infty.
\end{gathered}
\end{equation}
Using arguments similar to the proof of Proposition \ref{Syp11} we show that
 (a subsequence relabeled as) $(\hat{u}_{\delta n},\hat{v}_{\delta n})$
converges (on compact sets) to a solution $(\hat{u}_{\delta},\hat{v}_{\delta})$
of problem \eqref{NLmp4} with
$(\underline{u}_{\delta},\underline{v}_{\delta})\leq (\hat{u}_{\delta},
\hat{v}_{\delta}) \leq (\overline{u}_{\delta},\overline{v}_{\delta})$.
Likewise, (a subsequence of) $(\check{u}_{\delta n},\check{v}_{\delta n})$
converges (on compact sets) to a solution
$(\check{u}_{\delta},\check{v}_{\delta})$ of problem \eqref{NLmp4} with
$(\underline{u}_{\delta},\underline{v}_{\delta})
\leq (\check{u}_{\delta},\check{v}_{\delta}) \leq
(\overline{u}_{\delta},\overline{v}_{\delta})$.
Observe that, by \eqref{ineqf1} and \eqref{ineqf2},
$\check{f}_{1\delta}(x,t,\check{u}_{\delta}, \overline{v})
\geq \hat{f}_{1\delta}(x,t,\check{u}_{\delta}, \underline{v}) $,
$\check{f}_{2\delta}(x,t,\overline{u},\check{v}_{\delta})
\geq \hat{f}_{2\delta}(x,t,\underline{v}, \check{v}_{\delta}) $,
$\check{g}_{1\delta}(x,t,\check{u}_{\delta}, \overline{v})
\geq \hat{g}_{1\delta}(x,t,\check{u}_{\delta}, \underline{v}) $, and
$\check{g}_{2\delta}(x,t,\overline{u},\check{v}_{\delta})
\geq \hat{g}_{2\delta}(x,t,\underline{v}, \check{v}_{\delta}) $.
Therefore,
\begin{gather*}
\frac{\partial(\check{u}_{\delta}-\hat{u}_{\delta})}{\partial
t}-L(\check{u}_{\delta}-\hat{u}_{\delta})
\geq \hat{f}_{1\delta}(x,t,\check{u}_{\delta},
 \overline{v})-\hat{f}_{1\delta}(x,t,\hat{u}_{\delta}, \underline{v})\quad
\text{a.e. in } \Omega\times {\mathbb{R}},
\\
\mathcal{B}_\epsilon (\check{u}_{\delta}-\hat{u}_{\delta})
\geq \hat{g}_{1\delta}(x,t,\check{u}_{\delta}, \overline{v})
-\hat{g}_{1\delta}(x,t,\hat{u}_{\delta}, \underline{v})\quad \text{on }
\partial\Omega\times {\mathbb{R}},
\\
\frac{\partial(\check{v}_{\delta}-\hat{v}_{\delta})}{\partial
t}-L(\check{v}_{\delta}-\hat{v}_{\delta})
\geq \hat{f}_{2\delta}(x,t,\overline{u},\check{v}_{\delta})
- \hat{f}_{2,\delta}(x,t,\underline{v}, \hat{v}_{\delta}) \quad
\text{a.e. in } \Omega\times {\mathbb{R}},
\\
\mathcal{B}_\epsilon (\check{v}_{\delta}-\hat{v}_{\delta})
\geq \hat{g}_{2\delta}(x,t,\overline{u},\check{v}_{\delta})
- \hat{g}_{2\delta}(x,t,\underline{v}, \hat{v}_{\delta}) \quad\text{on }
\partial\Omega\times {\mathbb{R}}.
\end{gather*}
The monotonicity of the functions $\hat{f}_{i\delta}$, $\hat{g}_{i\delta}$,
$\check{f}_{i\delta}$  and $\check{g}_{i\delta} $, and an argument similar
to the one used in the proof of \cite[Proposition 2.7]{NM11}
imply that $\hat{u}_{\delta}\leq\check{u}_{\delta}$  and
$\hat{v}_{\delta}\leq\check{v}_{\delta}$ in
$\overline{\Omega}\times {\mathbb{R}}$. The proof is complete.
\end{proof}

\begin{lemma}\label{lemmodp}
Assume that {\rm (A1)--(A3)} are satisfied and that
$(\underline{u},\underline{v})$
and $(\overline{u},\overline{v})$ are subsolution and supersolution of
problem \eqref{Syeq1}
with $\underline{u}\leq \overline{u}$ and $\underline{v}\leq \overline{v}$.
Then there exist unique solutions
$(\hat{u},\hat{v}),  (\check{u},\check{v})\in
W^{2,1}_{p,\rm{loc}}(\Omega\times {\mathbb{R}})\cap
L^{\infty}({\Omega}\times {\mathbb{R}})\times W^{2,1}_{p,\rm{loc}}
(\Omega\times {\mathbb{R}})\cap
L^{\infty}({\Omega}\times {\mathbb{R}})$ to the respective problems
\eqref{NLmp2} and \eqref{NLmp1} such that
$$
\underline{u}\leq\hat{u}\le
\check{u} \leq\overline{u}\quad\text{and}\quad
\underline{v}\leq\hat{v}\le \check{v} \leq\overline{v}.
$$
\end{lemma}

\begin{proof}
Since the systems \eqref{NLmp2} and \eqref{NLmp1} are uncoupled, we have
that the uniqueness follows from the (nonincreasing) monotonicity
of the nonlinearities involved and Proposition \ref{eq1S}.
For $n\in \mathbb{N}$, let $\underline{u}_n=\underline{u}-\frac{1}{n}z$,
$\underline{v}_n=\underline{v}-\frac{1}{n}z$,
$\overline{u}_n=\overline{u}+\frac{1}{n}z$, and
$\overline{v}_n=\overline{v}+\frac{1}{n}z$, where $z$ is defined in
Lemma \ref{lemdelta}.
Consider the boundary value problems \eqref{NLmp4} and \eqref{NLmp3}
where the right hand sides are replaced by
\begin{gather*}
\hat{f}_{1n}(x,t, u,\underline{v}):= {f}_1(x,t,\underline{u},
\underline{v})+m_{1}(\underline{u}_{n}, u);\\
\hat{f}_{2n}(x,t, \underline{u}, {v}):= {f}_2(x,t,\underline{u},
\underline{v})+m_{2}(\underline{v}_{n}, v)
\\
\check{f}_{1n}(x,t, u,\overline{v}):={f}_1(x,t,\overline{u},
\overline{v})+m_{1}(\overline{u}_{n}, u);\\
\check{f}_{2n}(x,t, \overline{u},v):={f}_2(x,t,\overline{u},
\overline{v})+m_{2}(\overline{v}_{n}, v)
\\
\hat{g}_{1n}(x,t, u,\underline{v}):= g_1(x,t,\underline{u},
\underline{v})+\rho_{1}(\underline{u}_{n}- u); \\
\hat{g}_{2n}(x,t, \underline{u}, {v}):= g_2(x,t,\underline{u},
\underline{v})+\rho_{2}(\underline{v}_{n}- v)
\\
\check{g}_{1n}(x,t, u,\overline{v}):=g_1(x,t,\overline{u},
\overline{v})+\varrho_{1}(\overline{u}_{n}- u) ; \\
\check{g}_{2n}(x,t, \overline{u},v):=g_2(x,t,\overline{u}, \overline{v})
+\varrho_{2}(\overline{v}_{n}- v),
\end{gather*}
respectively.  By applying Lemma \ref{lemdelta} with
$\delta= \frac{1}{n}$, we get that, for each $n\in\mathbb{N}$ there exist
unique solutions $(\hat{u}_n,\hat{v}_n)$, $(\check{u}_n,\check{v}_n)\in
W^{2,1}_{p,\rm{loc}}(\Omega\times {\mathbb{R}})\cap
L^{\infty}({\Omega}\times {\mathbb{R}})\times W^{2,1}_{p,\rm{loc}}(\Omega\times {\mathbb{R}})\cap
L^{\infty}({\Omega}\times {\mathbb{R}})$ of these corresponding
problems such that
$$
\underline{u}_n\leq\hat{u}_n\le
\check{u}_n \leq\overline{u}_n\quad\text{and}\quad
\underline{v}_n\leq\hat{v}_n\le \check{v}_n \leq\overline{v}_n.
$$
Next, we show that the sequences $\{\hat{u}_n\}$, $\{\hat{v}_n\}$,
$\{\check{u}_n\}$,  and
$\{\check{v}_n\}$ are monotone and converge to unique solutions
of \eqref{NLmp2} and \eqref{NLmp1}, respectively. From the
definitions of $\underline{u}_n$, $\underline{v}_n$, $\overline{u}_n$,
and $\overline{v}_n$ we have that
the functions $\hat{f}_{1n}(x,t,\cdot,\underline{v})$ and
$\hat{g}_{1n}(x,t,\cdot,\underline{v})$ are nonincreasing for
$u\in[\underline{u}_n,\infty)$, and $\hat{f}_{2n}(x,t,\underline{u},\cdot)$
 and $\hat{g}_{2n}(x,t,\underline{u},\cdot)$ are nonincreasing for
$v\in[\underline{v}_n,\infty)$, and that they are nondecreasing with
respect to $n$. Similarly,  the functions
$\check{f}_{1n}(x,t,\cdot,\overline{v})$ are nonincreasing for
$u\in(-\infty,\overline{v}_n]$, and
$\check{g}_{1n}(x,t,\cdot,\overline{v})$ are nonincreasing for
$\check{f}_{2n}(x,t,\overline{u},\cdot)$ and
$\check{g}_{2n}(x,t,\overline{u},\cdot)$ are nonincreasing for
$v\in(-\infty,\overline{v}_n]$,  and they are nonincreasing with
respect to $n$.
Therefore,
\begin{gather*}
\frac{\partial(\hat{u}_{n}-\hat{u}_{n-1})}{\partial t}-L(\hat{u}_{n}
-\hat{u}_{n-1}) \geq \hat{f}_{1(n-1)}(x,t,\hat{u}_{n}, \underline{v})
-\hat{f}_{1(n-1)}(x,t,\hat{u}_{n-1}, \underline{v}),
\\
\mathcal{B}_\epsilon (\hat{u}_{n}-\hat{u}_{n-1})
\geq \hat{g}_{1(n-1)}(x,t,\hat{u}_{n}, \underline{v})
-\hat{g}_{1(n-1)}(x,t,\hat{u}_{n-1}, \underline{v}),
\\
\frac{\partial(\hat{v}_{n}-\hat{v}_{n-1})}{\partial
t}-L(\hat{v}_{n}-\hat{v}_{n-1})
\geq \hat{f}_{2(n-1)}(x,t,\underline{u} ,\hat{v}_{n})
 -\hat{f}_{2(n-1)}(x,t,\underline{u}, \hat{v}_{n-1}),
\\
\mathcal{B}_\epsilon (\hat{v}_{n}-\hat{v}_{n-1})
\geq \hat{g}_{2(n-1)}(x,t,\underline{u} ,\hat{v}_{n})
 -\hat{g}_{2(n-1)}(x,t,\underline{u}, \hat{v}_{n-1}).
\end{gather*}
By the monotonicity of $\hat{f}_{1n}(x,t,\cdot,\underline{v}),
\hat{g}_{1n}(x,t,\cdot,\underline{v})$,
$\hat{f}_{2n}(x,t,\underline{u},\cdot)$,
$\hat{g}_{2n}(x,t,\underline{u},\cdot)$  and Corollary \ref{Sycor2},
we show as in the proof of Lemma \ref{lemdelta} that
$\hat{u}_{n-1}\leq\hat{u}_{n}$ and
$\hat{v}_{n-1}\leq\hat{v}_{n}$. Similarly, we get that
$\check{u}_{n}\leq\check{u}_{n-1}$ and $\check{v}_{n}\leq\check{v}_{n-1}$.
Therefore,
\begin{gather*}
\underline{u}_1\leq\hat{u}_1\leq \hat{u}_2\leq \dots\leq
\hat{u}_{n-1}\leq \hat{u}_n \leq \ldots\leq
\check{u}_n\leq \check{u}_{n-1}\leq \ldots\leq
\check{u}_2\leq \check{u}_1 \leq \overline{u}_1.
\\
\underline{v}_1\leq\hat{v}_1\leq \hat{v}_2\leq \dots\leq
\hat{v}_{n-1}\leq \hat{v}_n \leq \ldots\leq
\check{v}_n\leq \check{v}_{n-1}\leq \ldots\leq
\check{v}_2\leq \check{v}_1 \leq \overline{v}_1.
\end{gather*}
It follows that $\{(\hat{u}_n,\hat{v}_n)\}$ and
$\{(\check{u}_n,\check{v}_n)\}$ converge
(pointwise) to $(\hat{u},\hat{v})$ and $(\check{u},\check{v})$, respectively,
with $(\underline{u},\underline{v})\leq (\hat{u},\hat{v})\leq (\check{u},\check{v})\leq
(\overline{u},\overline{v})$.
Now, we will show that $\underline{v}$ and $\overline{v}$ are
respective solutions of \eqref{NLmp2} and \eqref{NLmp1}.
Consider
$Q_1=\Omega\times (-1, 1)$ and $Q_2= \Omega\times (-2, 2). $ For each
$n\in {\mathbb{N}}$, define $\hat{z}_n(x,t)= \zeta(t) \hat{u}_{ n}(x,t)$,
$\hat{w}_n(x,t)= \zeta(t) \hat{v}_{n}(x,t)$,
for all $(x, t)\in \overline{\Omega}\times [-2,  2]$, where
$\zeta\in C^{\infty}({\mathbb{R}}),  0\leq \zeta\leq 1$ and
$\zeta(s)=0$ if $s\leq -2,  \zeta(s)=1\;
\text{if}\;s\geq -(2-\delta)$ with $0<\delta<1$. Observe that
$\hat z_n=\hat u_n $  and $\hat w_n= \hat v_n$, in
$\overline{\Omega}\times [-1,  1]$,
and satisfy the  linear uncoupled system
\begin{equation}\label{Syeq4c}
\begin{gathered}
 \frac{\partial {\hat{z}}_n}{\partial t}-L_1{\hat{z}}_n
= \frac {d\zeta}{dt}\hat{u}_{ n}+\zeta \hat{f}_{1n}(x,t,\hat u_n,\underline{v})
\quad\text{in } \Omega\times (-2, 2],\\
\frac{\partial {\hat{w}}_n}{\partial t}-L_2{\hat{w}}_n
=\frac {d\zeta}{dt} \hat{v}_{ n} +\zeta \hat{f}_{2n}(x,t,\overline{u},\hat v_n)
 \quad\text{in } \Omega\times (-2, 2],\\
\mathcal{B}_1 {\hat{z}}_n = \zeta \hat{g}_{1n}(x,t,\hat u_n,\underline{v})
 \quad\text{on } \partial \Omega\times(-2, 2],\\
\mathcal{B}_2 {\hat{w}}_n =\zeta \hat{g}_{2n}(x,t,\overline{u},\hat v_n)
 \quad\text{on }\partial \Omega\times(-2, 2],
\\
\hat{z}_n(x,-2)=0 \quad \text{in }\overline{\Omega},\\
\hat{ w}_n(x,-2)=0 \quad \text{in } \overline{\Omega},\\
 \sup_{\Omega\times {\mathbb{R}}}\{|{\hat{z}}_n(x,t)|,
|{\hat{w}}_n(x,t)|\}<\infty.
\end{gathered}
\end{equation}
Arguments similar to those used in the proof of Proposition \ref{Syp11}
show that   $(\hat{u},\hat{v})$ and
$(\check{u},\check{v})\in W_{p,{\rm loc}}^{2,1}({\Omega}\times
{\mathbb{R}})\cap L^{\infty}({\Omega}\times {\mathbb{R}})
\times W_{p,{\rm loc}}^{2,1}({\Omega}\times
{\mathbb{R}})\cap L^{\infty}({\Omega}\times {\mathbb{R}})$  are
solutions of \eqref{NLmp2} and \eqref{NLmp1}, respectively,
with $(\underline{u},\underline{v})\leq (\hat{u},\hat{v})
\leq (\check{u},\check{v})\leq (\overline{u},\overline{v})$.
 The proof is complete.
\end{proof}

\begin{proof}[Proof of Theorem \ref{Syt1}]
We construct two sequences $\{(\hat{u}_n,\hat{v}_n)\}$ and
$\{(\check{u}_n,\check{v}_n)\}$  successively from the
nonlinear iteration processes

\begin{equation}\label{mhp1}
\begin{gathered}
 \frac{\partial \hat{u}_n}{\partial t}(x,t)-L\hat{u}_n (x,t)=
 f_1(x,t,\hat{u}_{n-1},\hat{v}_{n-1})+m_1(\hat{u}_{n-1},\hat{u}_n)
\quad\text{a.e. in } \Omega\times {\mathbb{R}},\\
 \frac{\partial \hat{v}_n}{\partial t}(x,t)-L\hat{v}_n (x,t)
=  f_2(x,t,\hat{u}_{n-1},\hat{v}_{n-1})+m_2(\hat{v}_{n-1},\hat{v}_n)
 \quad\text{a.e. in } \Omega\times {\mathbb{R}},\\
\mathcal{B}_1 \hat{u}_n = g_1(x,t,\hat{u}_{n-1},\hat{v}_{n-1})
+\varrho_1(\hat{u}_{n-1}-\hat{u}_n)
\quad \text{on } \partial \Omega\times {\mathbb{R}},\\
\mathcal{B}_2 \hat{v}_n =
g_2(x,t,\hat{u}_{n-1},\hat{v}_{n-1})+\varrho_2(\hat{v}_{n-1}-\hat{v}_n)
\quad \text{on } \partial \Omega\times {\mathbb{R}},\\
 \sup_{\Omega\times
{\mathbb{R}}}\{|\hat{u}_n(x,t)|,|\hat{v}_n(x,t)|\}<\infty,
\end{gathered}
\end{equation}
and
\begin{equation}\label{mhp2}
\begin{gathered}
 \frac{\partial \check{u}_n}{\partial t}(x,t)-L\check{u}_n (x,t)
= f_1(x,t,\check{u}_{n-1},\check{v}_{n-1})+m_1(\check{u}_{n-1},\check{u}_n)
 \quad\text{a.e. in } \Omega\times {\mathbb{R}},\\
 \frac{\partial \check{v}_n}{\partial t}(x,t)-L\check{v}_n (x,t)=
 f_2(x,t,\check{u}_{n-1},\check{v}_{n-1})+m_2(\check{v}_{n-1},\check{v}_n)
 \quad\text{a.e. in } \Omega\times {\mathbb{R}},\\
\mathcal{B}_1 \check{u}_n =
g_1(x,t,\check{u}_{n-1},\check{v}_{n-1})+\varrho_1(\check{u}_{n-1}-\check{u}_n)
\quad \text{on } \partial \Omega\times {\mathbb{R}},\\
\mathcal{B}_2 \check{v}_n =
g_2(x,t,\check{u}_{n-1},\check{v}_{n-1})+\varrho_2(\check{v}_{n-1}-\check{v}_n)
\quad \text{on } \partial \Omega\times {\mathbb{R}},\\
\sup_{\Omega\times {\mathbb{R}}}\{|\check{u}_n(x,t)|,
|\check{v}_n(x,t)|\}<\infty.
\end{gathered}
\end{equation}
We show that these sequences are well defined and that they converge
monotonically to a solution of \eqref{Syeq1}.
Indeed, set
\begin{gather*} %\label{hfnlbar}
\hat{f}_{1n}(x,t,u,\hat{v}_{n-1})=f_1(x,t,\hat{u}_{n-1},\hat{v}_{n-1})
+m_1(\hat{u}_{n-1},u),
\\
 \hat{f}_{2n}(x,t,\hat{u}_{n-1},v)=f_2(x,t,\hat{u}_{n-1},\hat{v}_{n-1})
+m_2(\hat{v}_{n-1},v),
\\ %\label{hfnlbar}
\hat{g}_{1n}(x,t,u,\hat{v}_{n-1})=f_1(x,t,\hat{u}_{n-1},\hat{v}_{n-1})
+m_1(\hat{u}_{n-1},u),
\\
 \hat{g}_{2n}(x,t,\hat{u}_{n-1},v)=f_2(x,t,\hat{u}_{n-1},\hat{v}_{n-1})
+m_2(\hat{v}_{n-1},v),
\\ %\label{hfnlbar}
\check{f}_{1n}(x,t,u,\check{v}_{n-1})=f_1(x,t,\check{u}_{n-1},\check{v}_{n-1})
+m_1(\check{u}_{n-1},u),
\\
 \check{f}_{2n}(x,t,\check{u}_{n-1},v)=f_2(x,t,\check{u}_{n-1},\check{v}_{n-1})
+m_2(\check{v}_{n-1},v),
\\ %\label{hfnlbar}
\check{g}_{1n}(x,t,u,\check{v}_{n-1})=f_1(x,t,\check{u}_{n-1},\check{v}_{n-1})
+m_1(\check{u}_{n-1},u),
\\
 \check{g}_{2n}(x,t,\check{u}_{n-1},v)=f_2(x,t,\check{u}_{n-1},\check{v}_{n-1})
+m_2(\check{v}_{n-1},v),
\end{gather*}
where $(\hat{u}_0,\hat{v}_0)=(\underline{u},\underline{v})$ and
$(\check{u}_0,\check{v}_0)=(\overline{u},\overline{v})$.
It follows immediately from Lemma \ref{lemmodp} that the first
iterations $(\hat{u}_1,\hat{v}_1)$ in \eqref{mhp1} and
$(\check{u}_1,\check{v}_1)$ in \eqref{mhp2} exist and satisfy the inequalities
$(\underline{u},\underline{v})\leq(\hat{u}_1,\hat{v}_1)
\leq (\check{u}_1,\check{v})\leq (\overline{u},\overline{v})$,
when one starts with $(\hat{u}_0,\hat{v}_0)=(\underline{u},\underline{v})$
and $(\check{u}_0,\check{v}_0)=(\overline{u},\overline{v})$. For $n\geq 2$,
we use an induction argument to show that
$(\hat{u}_{n-1},\hat{v}_{n-1})\leq(\hat{u}_{n},\hat{v}_{n})\leq
(\check{u}_{n-1},\check{v}_{n-1})\leq (\check{u}_{n},\check{v}_{n})$.
However, in order to apply Lemma \ref{lemmodp} inductively, we need to
show that, at each iteration, the functions $(\hat{u}_{n-1},\hat{v}_{n-1})$
and $(\check{u}_{n-1},\check{v}_{n-1})$ are ordered subsolution and
supersolution of problem \eqref{Syeq1}.
By using \eqref{transfo2}, (A2), quasimonotonicity decreasing and the equations
\eqref{mhp1} and \eqref{mhp2}, we get
\begin{gather*}
\begin{aligned}
&\frac{\partial \hat{u}_{n-1}}{\partial t}-L\hat{u}_{n-1}
- f_1(x,t,\hat{u}_{n-1},\hat{v}_{n-1})\\
&= f_1(x,t,\hat{u}_{n-2},\hat{v}_{n-2})+m_1(\hat{u}_{n-2},\hat{u}_{n-1})
 -f_1(x,t,\hat{u}_{n-1},\hat{v}_{n-1})\leq 0,
\end{aligned} \\
\begin{aligned}
&\frac{\partial \hat{v}_{n-1}}{\partial t}-L\hat{v}_{n-1}
- f_2(x,t,\hat{u}_{n-1},\hat{v}_{n-1})\\
&= f_2(x,t,\hat{u}_{n-2},\hat{v}_{n-2})+m_2(\hat{v}_{n-2},\hat{v}_{n-1})
 -f_2(x,t,\hat{u}_{n-1},\hat{v}_{n-1}) \leq 0,
\end{aligned}\\
\begin{aligned}
&\mathcal{B} \hat{u}_{n-1}-g_1(x,t,\hat{u}_{n-1},\hat{v}_{n-1})\\
& = g_1(x,t,\hat{u}_{n-2},\hat{v}_{n-2})+ \varrho
(\hat{u}_{n-2}-\hat{u}_{n-1})
 -g_1(x,t,\hat{u}_{n-1},\hat{v}_{n-1}) \leq 0,
\end{aligned}\\
\begin{aligned}
&\mathcal{B} \hat{v}_{n-1}-g_2(x,t,\hat{u}_{n-1},\hat{v}_{n-1})\\
&= g_2(x,t,\hat{u}_{n-2},\hat{v}_{n-2})+ \varrho(\hat{v}_{n-2}-\hat{v}_{n-1})
-g_2(x,t,\hat{u}_{n-1},\hat{v}_{n-1}) \leq 0,
\end{aligned}\\
 \sup_{\Omega\times
{\mathbb{R}}}\{|\hat{u}_{n-1}(x,t)|,|\hat{v}_{n-1}(x,t)|\}<\infty,
\end{gather*}
 and
\begin{gather*}
\begin{aligned}
&\frac{\partial \check{u}_{n-1}}{\partial t}-L\check{u}_{n-1}
- f_1(x,t,\check{u}_{n-1},\check{v}_{n-1})\\
&= f_1(x,t,\check{u}_{n-2},\check{v}_{n-2})+m_1(\check{u}_{n-2},\check{u}_{n-1})
  -f_1(x,t,\check{u}_{n-1},\check{v}_{n-1})\geq 0,
\end{aligned}\\
\begin{aligned}
& \frac{\partial \check{v}_{n-1}}{\partial t}-L\check{v}_{n-1}
 - f_2(x,t,\check{u}_{n-1},\check{v}_{n-1})\\
&=  f_2(x,t,\check{u}_{n-2},\check{v}_{n-2})+m_2(\hat{v}_{n-2},\hat{v}_{n-1})
 -f_2(x,t,\check{u}_{n-1},\check{v}_{n-1})\geq 0,
\end{aligned}\\
\begin{aligned}
&\mathcal{B} \check{u}_{n-1}-g_1(x,t,\check{u}_{n-1},\check{v}_{n-1})\\
& = g_1(x,t,\check{u}_{n-2},\check{v}_{n-2})+ \varrho_1
(\check{u}_{n-2}-\check{u}_{n-1})
 -g_1(x,t,\check{u}_{n-1},\check{v}_{n-1}) \geq 0,
\end{aligned}\\
\begin{aligned}
&\mathcal{B} \check{v}_{n-1}-g_2(x,t,\check{u}_{n-1},\check{v}_{n-1})\\
&=  g_2(x,t,\check{u}_{n-2},\check{v}_{n-2})+ \varrho_2
(\check{v}_{n-2}-\check{v}_{n-1})
-g_2(x,t,\check{u}_{n-1},\check{v}_{n-1}) \geq 0,
\end{aligned}\\
 \sup_{\Omega\times
{\mathbb{R}}}\{|\check{u}_{n-1}(x,t)|,|\check{v}_{n-1}(x,t)|\}<\infty,
\end{gather*}
which shows that,  for $n\geq 2$, the functions
$(\hat{u}_{n-1},\hat{v}_{n-1})$ and $(\check{u}_{n-1},\check{v}_{n-1})$
are ordered subsolution and supersolution of \eqref{Syeq1}.
Since
$\hat{f}_{1n},\hat{f}_{2n}$, $\check{f}_{1n},\check{f}_{2n}$, and
$\hat{g}_{1n},\hat{g}_{2n}$, $\check{g}_{1n},\check{g}_{2n}$
satisfy the conditions of Lemma \ref{lemmodp} with
$(\hat{u}_{n-1},\hat{v}_{n-1})$ as subsolution and
$(\check{u}_{n-1},\check{v}_{n-1})$ as
supersolution, the existence of solutions $(\hat{u}_{n},\hat{v}_{n})$ to
Eq.\eqref{mhp1} and $(\check{u}_{n},\check{v}_{n})$ to Eq.\eqref{mhp2}
such that $(\hat{u}_{n-1},\hat{v}_{n-1})\leq(\hat{u}_{n},\hat{v}_{n})
\leq (\check{u}_{n},\check{v}_{n})\leq (\check{u}_{n-1},\check{v}_{n-1})$
is ensured by Lemma \ref{lemmodp}.
We therefore have that
\begin{gather*}
\underline{u}_1\leq\hat{u}_1\leq \hat{u}_2\leq \dots\leq
\hat{u}_{n-1}\leq \hat{u}_n \leq \ldots\leq
\check{u}_n\leq \check{u}_{n-1}\leq \ldots\leq
\check{u}_2\leq \check{u}_1 \leq \overline{u}_1,
\\
\underline{v}_1\leq\hat{v}_1\leq \hat{v}_2\leq \dots\leq
\hat{v}_{n-1}\leq \hat{v}_n \leq \ldots\leq
\check{v}_n\leq \check{v}_{n-1}\leq \ldots\leq
\check{v}_2\leq \check{v}_1 \leq \overline{v}_1.
\end{gather*}
It follows that $\{(\hat{u}_n,\hat{v}_n)\}$ and $\{(\check{u}_n,\check{v}_n)\}$
converge (pointwise) to $(\hat{u},\hat{v})$ and $(\check{u},\check{v})$,
respectively, with
$(\underline{u},\underline{v})\leq (\hat{u},\hat{v})
\leq (\check{u},\check{v})\leq (\overline{u},\overline{v})$.
Consider  $Q_1=\Omega\times (-1, 1)$ and $Q_2= \Omega\times (-2, 2)$.
 For each $n\in
{\mathbb{N}}$, define $\hat{z}_n(x,t)= \zeta(t) \hat{u}_{ n}(x,t)$,
 $\hat{w}_n(x,t)= \zeta(t) \hat{v}_{n}(x,t)$,
for all $(x, t)\in \overline{\Omega}\times [-2,  2]$, where
$\zeta\in C^{\infty}({\mathbb{R}})$, $ 0\leq\zeta\leq 1$ and
 $\zeta(s)=0$  \text{if } $s\leq -2$,  $\zeta(s)=1$
if $s\geq -(2-\delta)$ with $0<\delta<1$. Observe that
$\hat z_n=\hat u_n $  and $\hat w_n= \hat v_n$, in
$\overline{\Omega}\times [-1,  1]$ and $(z_n,w_n)$ satisfies the following
uncouple system
\begin{gather*}
 \frac{\partial {\hat{z}}_n}{\partial t}-L_1{\hat{z}}_n = \frac
{d\zeta}{dt}\hat{u}_{ n}+\zeta \hat{f}_{1n}(x,t,\hat u_n,\hat{v}_{n-1})
\quad\text{in } \Omega\times (-2, 2],\\
\frac{\partial {\hat{w}}_n}{\partial t}-L_2{\hat{w}}_n =\frac
{d\zeta}{dt} \hat{v}_{ n} +\zeta \hat{f}_{2n}(x,t,\hat{u}_{n-1},\hat v_n)
\quad\text{in } \Omega\times (-2, 2],\\
\mathcal{B}_1 {\hat{z}}_n = \zeta \hat{g}_{1n}(x,t,\hat u_n,\hat{v}_{n-1})
\quad\text{on } \partial \Omega\times (-2, 2],\\
\mathcal{B}_2 {\hat{w}}_n =\zeta \hat{g}_{2n}(x,t,\hat u_{n-1},\hat{v}_{n})
\quad\text{on } \partial \Omega\times (-2, 2],\\
\hat{z}_n(x,-2)=0 \quad \text{in }\overline{\Omega},\\
\hat{ w}_n(x,-2)=0 \quad \text{in }\overline{\Omega},\\
 \sup_{\Omega\times {\mathbb{R}}}\{|{\hat{z}}_n(x,t)|,
|{\hat{w}}_n(x,t)|\}<\infty,
\end{gather*}
where
$(z_n,\, w_n)\in
\mathrm{W}_p^{2,  1}(Q_2)\times \mathrm{W}_p^{2,  1}(Q_2) $
(with $p=\frac{N+2}{1-\mu}$). Moreover
\begin{gather*}
|z_n|_{W_p^{2,  1}(Q_2)}\leq K_0\Big(|\frac
{d\zeta}{dt}\hat{u}_n+\zeta \hat{f}_{1,n}|_{L^p(Q_2)}+|\zeta
\hat{g}_{1,n}|_{W_p^{2-\epsilon-\frac{1}{p}, (2-\epsilon-\frac{1}{p})/2}({\partial
\Omega\times(-2, 2)})}\Big),
\\
|w_n|_{W_p^{2,  1}(Q_2)}\leq K_0\Big(|\frac
{d\zeta}{dt}\hat{v}_n+\zeta \hat{f}_{2,n}|_{L^p(Q_2)}+|\zeta
\hat{g}_{2n}|_{W_p^{2-\epsilon-\frac{1}{p}, (2-\epsilon-\frac{1}{p})/2}({\partial
\Omega\times(-2, 2)})}\Big),
 \end{gather*}
for all $n\in {\mathbb{N}}$,  where $K_0$ is a constant which depends on
$Q_2$. Set $V_n=(z_n,w_n)$  with $|V_n|_{W_p^{2,
1}(Q_2)}=|z_n|_{W_p^{2,  1}(Q_2)}+|w_n|_{W_p^{2,  1}(Q_2)}$.
Observe that for the Dirichlet boundary condition, we get immediately that
$|V_n|_{W_p^{2,  1}(Q_2)}\leq { C}$, for all $n$. To show that
$|V_n|_{W_p^{2,  1}(Q_2)}\leq {C}$  for all $n$ for the Neumann boundary
condition,  we proceed as follows. Using assumptions (A2) we compute
$|\zeta\hat{g}_{i\,n}|_{W_p^{1-\frac{1}{p}, (1-\frac{1}{p})/2}({\partial
\Omega\times(-2, 2)})}$ ($i=1,2$) to obtain
\begin{align*} %\label{SylPsi3}
&|\zeta \hat{g}_{i\,n}|_{W_p^{1-\frac{1}{p}, (1-\frac{1}{p})/2}({\partial
\Omega\times(-2, 2)})} \\
&\leq \hat{C} (1+|V_{n}|_{W_p^{1-\frac{1}{p}, (1-\frac{1}{p})/2}({\partial
\Omega\times(-2, 2)})}
+|V_{n-1}|_{W_p^{1-\frac{1}{p}, (1-\frac{1}{p})/2}({\partial
\Omega\times(-2, 2)})}),
\end{align*}
where $\hat{C}$ is independent of $n$ since
$|\zeta \hat{g}_{i,n}|_{L^p({\partial \Omega\times(-2, 2)})}\leq \text{const}$
for all $n \in {\mathbb{N}}$.
Using the continuity of the trace operator and the fact that
$\hat{u}_n, \hat{v}_n$ and $\hat{f}_{i,n}$ are (uniformly) bounded,
we get that
$$
|V_n|_{W_p^{2,  1}(Q_2)}\leq \tilde{K}
\Big(1+|V_{n}|_{W_p^{1,1/2}(Q_2)}+|V_{n-1}|_{W_p^{{1,1/2}(Q_2)}}\Big),
$$
where ${K}$ is independent of $n$.
Using the interpolation inequality  \eqref{Syintp2}, we obtain
\begin{equation*} %\label{SylPsi5}
|V_n|_{W_p^{2,  1}(Q_2)}\leq \frac{K}{1-\varepsilon K}
\Big(1+\frac{1}{\varepsilon}|V_{n}|_{L^p(Q_2)}+
|V_{n-1}|_{W_p^{1, 1/2}(Q_2)}\Big),
\end{equation*}
where $\varepsilon\tilde{K}<1$. Since ${V}_n$ is
(uniformly) bounded, we deduce that
\begin{equation*} %\label{SylPsi5}
|V_n|_{W_p^{2,  1}(Q_2)}\leq C \Big(1+|V_{n-1}|_{W_p^{1, 1/2}(Q_2)}\Big),
\end{equation*}
where $C$ is a constant independent of $n$.
Using the same reasoning as in the proof of Proposition \ref{Syp11},
one shows that $\{V_n\}$ converges to a solution $(\hat{u},\hat{v})$
of \eqref{Syeq1} with $(\hat{u},\hat{v})\in W_{p,{\rm loc}}^{2,1}(\Omega\times
{\mathbb{R}})\cap L^{\infty}({\Omega}\times {\mathbb{R}})\times W_{p,{\rm loc}}^{2,1}(\Omega\times
{\mathbb{R}})\cap L^{\infty}({\Omega}\times {\mathbb{R}}) $.  An
analogous argument shows also that
$(\check{u},\check{v})\in W_{p,{\rm loc}}^{2,1}(\Omega\times
{\mathbb{R}})\cap L^{\infty}({\Omega}\times {\mathbb{R}})
\times W_{p,{\rm loc}}^{2,1}(\Omega\times
{\mathbb{R}})\cap L^{\infty}({\Omega}\times {\mathbb{R}}) $
is a solution of \eqref{Syeq1}, and  that
$(\underline{u},\underline{v})\leq (\hat{u},\hat{v})
\leq (\check{u},\check{v})\leq
(\overline{u},\overline{v})$.
\end{proof}

\section{Examples}

In this section, we illustrate our results with the following examples.

\begin{example}\rm (Cooperative Model)
Consider the system
\begin{equation}{\label{Syexam1}}
\begin{gathered}
 \frac{\partial u}{\partial t}-\triangle u =
u(a_1(x,t)-b_1(x,t) u+c_1(x,t) v) \quad\text{in } \Omega\times
{\mathbb{R}},\\
\frac{\partial v}{\partial t}-\triangle v = v(a_2(x,t)+b_2(x,t)
u-c_2(x,t) v)\quad\text{in } \Omega\times
{\mathbb{R}},\\
u=0=v \quad\text{on } \partial{\Omega}\times {\mathbb{R}},\\
 \sup_{\Omega\times
{\mathbb{R}}}\{|u|,|v|\}<\infty,
\end{gathered}
\end{equation}
where $ a_i:\Omega\times {\mathbb{R}}\to {\mathbb{R}}$,
$b_i:\Omega\times {\mathbb{R}}\to {\mathbb{R}}$,
$c_i:\Omega\times {\mathbb{R}}\to {\mathbb{R}}$  are in
$C_{\rm{loc}}^{\mu,\mu/2}(\overline{\Omega}\times {\mathbb{R}})$. 
The  coefficients satisfy the following conditions:
For all $(x,t)\in\Omega\times {\mathbb{R}}$,  
$\lambda_1 <\alpha_i\leq a_i(x,t)\leq A_i$,
$0<\beta_i\leq b_i(x,t)\leq B_i$,
$$ 
0<\gamma_i\leq c_i(x,t)\leq C_i, \quad\text{and} \quad
\frac{B_2}{\beta_1}<\frac{\gamma_2}{C_1}.
$$
Here, $\lambda_1$ is the first eigenvalue of the Laplacian, and
$\alpha_i, \beta_i, \gamma_i, A_i, B_i, C_i \in {\mathbb{R}}$.
Observe that the presence of the $u-$population species encourages
the growth of the $v-$population species and vice versa. So, the
reaction functions $f_1(x,t,u,v)= u(a_1(x,t)-b_1(x,t) u+c_1(x,t) v)$
and $f_2(x,t,u,v)= v(a_2(x,t)+b_2(x,t) u-c_2(x,t) v)$ are
quasimonotone nondecreasing in $[0,\infty)\times[0,\infty)$. In
order to apply Theorem {\ref{Syt1}}, we need to have an ordered
sub-solution $(\underline{u},\underline{v})$ and  super-solution
$(\overline{u},\overline{v})$ of \eqref{Syexam1} that  satisfy the
following inequalities
\begin{gather*} %{\label{coop1}}
\frac{\partial \underline{u}}{\partial t}-\triangle
\underline{u} \leq \underline{u}(a_1(x,t)-b_1(x,t)
\underline{u}+c_1(x,t) \underline{v})  \quad\text{in } \Omega\times
{\mathbb{R}},\\
\frac{\partial \underline{v}}{\partial t}-\triangle \underline{v}
\leq \underline{v}(a_2(x,t)+b_2(x,t) \underline{u}-c_2(x,t)
\underline{v})\quad\text{in } \Omega\times
{\mathbb{R}},\\
 \frac{\partial \overline{u}}{\partial t}-\triangle
\overline{u} \ge \overline{u}(a_1(x,t)-b_1(x,t)
\overline{u}+c_1(x,t) \overline{v})  \quad\text{in } \Omega\times
{\mathbb{R}},\\
 \frac{\partial \overline{v}}{\partial
t}-\triangle \overline{v} \geq \overline{v}(a_2(x,t)+b_2(x,t)
\overline{u}-c_2(x,t) \overline{v}  \quad\text{in } \Omega\times
{\mathbb{R}},\\
\underline{u}\le 0\le \overline{u} \quad\text{on } 
\partial{\Omega}\times {\mathbb{R}},\\
\underline{v} \leq 0\leq  \overline{v} \quad\text{on } 
\partial{\Omega}\times {\mathbb{R}},\\
 {{\sup_{\Omega\times
{\mathbb{R}}}\{|\underline{u}|,|\overline{u}|,|\underline{v}|,|\overline{v}|\}
<\infty.}}
\end{gather*}
  Choose $(\underline{u},\underline{v})= (\varepsilon
\varphi, \varepsilon \varphi)$
 with $0<\varepsilon<\min\{\frac{\alpha_1 -\lambda_1}{B_1}, 
\frac{\alpha_2 -\lambda_1}{B_2}\}$
 where $0<\varphi$ is the eigenfunction associated with $\lambda_1$. 
Pick $(\overline{u},\overline{v})= (M_1,M_2)$, where 
$ M_1=\tau\frac{A_1\gamma_2+A_2C_1}{\beta_1\gamma_2-B_2C_1}$
 and $ M_2=\tau\frac{A_2\beta_1+A_1B_2}{\beta_1\gamma_2-B_2C_1}$. The
  positive constant
 $\tau$ is  chosen so that $(\varepsilon \varphi,
\varepsilon \varphi)\leq (M_1, M_2)$. Then by Theorem {\ref{Syt1}},
there exists  a positive solution $(u,v)$ such that 
$(\varepsilon \varphi, \varepsilon \varphi)\leq (u,v)\leq (M_1, M_2)$
in $\overline{\Omega}\times {\mathbb{R}}$. Thus, $(u,v)$ does
not tend to zero as $t$ tends to $\pm \infty$, for each
$x\in\Omega$.
\end{example}

\begin{example} \rm (Generalized Cooperative Model with nonlinear Boundary
 conditions)
\begin{equation}\label{Syexam2}
\begin{gathered}
\frac{\partial u}{\partial t}-\triangle u 
= u^m(a_1(x,t)-b_1(x,t) u+c_1(x,t) v) \quad\text{in } \Omega\times
{\mathbb{R}},\\
\frac{\partial v}{\partial t}-\triangle v = v^m(a_2(x,t)+b_2(x,t)
u-c_2(x,t) v)\quad\text{in } \Omega\times
{\mathbb{R}},\\
\frac{\partial u}{\partial\nu} =u^n\,(\delta_1-u +\sigma_1 v) \quad\text{on } \partial{\Omega}\times
{\mathbb{R}},\\
\frac{\partial v}{\partial\nu} =v^n\,(\delta_2-v +\sigma_2 u) \quad\text{on }
 \partial{\Omega}\times {\mathbb{R}},\\
 \sup_{\Omega\times
{\mathbb{R}}}\{|u|,|v|\}<\infty,
\end{gathered}
\end{equation}
where $ n,m\in {\mathbb{N}},\,0<\delta_i\in {\mathbb{R}}$,
$0<\sigma_i<1$,
\begin{equation*} 
a_i:\Omega\times
{\mathbb{R}}\to {\mathbb{R}},b_i:\Omega\times
{\mathbb{R}}\to {\mathbb{R}}, c_i:\Omega\times
{\mathbb{R}}\to {\mathbb{R}}\quad\text{are in }
C_{\rm{loc}}^{\mu,\,\mu/ 2}(\overline{\Omega}\times {\mathbb{R}}).
\end{equation*}
For all $(x,t)\in\Omega\times {\mathbb{R}}$, 
$ 0 <\alpha_i\leq a_i(x,t)\leq A_i$, 
$ 0<\beta_i\leq b_i(x,t)\leq B_i$,
$ 0<\gamma_i\leq c_i(x,t)\leq C_i$,  where 
$\alpha_i, \beta_i, \gamma_i,  A_i, B_i, C_i \in {\mathbb{R}}$.
 Note that the nonlinearities satisfy the quasimonotone
nondecreasing property.
Choose $(\underline{u},\underline{v})= (D, D)$, where $D>0$ is very
small such that $D<\min\{\delta_1, \delta_2\}$. Pick
$(\overline{u},\overline{v})= (M,M)$, where $M$ is a constant such
that $M>\max\left\{\frac{\delta_1}{1-\sigma_1},
\frac{\delta_2}{1-\sigma_2}\right\}$.\\
Then it follows from Theorem \ref{Syt1} that the system
\eqref{Syexam2} has a positive solution $(u,v)$ such that 
$(D, D)\leq (u,v)\leq (M, M)$ in $\overline{\Omega}\times
{\mathbb{R}}$. Thus, $(u,v)$ does not tend to zero as $t$ tends to
$\pm \infty$.
\end{example}

\begin{example}\rm (Competitive Model)
\begin{equation}\label{Syexam3}
\begin{gathered}
 \frac{\partial u}{\partial t}-\triangle u =
u(a_1(x,t)-b_1(x,t) u-c_1(x,t) v) \quad\text{in } \Omega\times
{\mathbb{R}},\\
\frac{\partial v}{\partial t}-\triangle v = v(a_2(x,t)-b_2(x,t)
v-c_2(x,t) u)\quad\text{in } \Omega\times
{\mathbb{R}},\\
\frac{\partial u}{\partial\nu}=0=\frac{\partial v}{\partial\nu} 
\quad\text{on } \partial{\Omega}\times {\mathbb{R}},\\
 \sup_{\Omega\times
{\mathbb{R}}}\{|u|,|v|\}<\infty,
\end{gathered}
\end{equation}
where $a_i:\Omega\times{\mathbb{R}}\to {\mathbb{R}}$,
$b_i:\Omega\times {\mathbb{R}}\to {\mathbb{R}}$, 
$c_i:\Omega\times {\mathbb{R}}\to {\mathbb{R}}$  are in 
$C_{\rm{loc}}^{\mu,\,\mu/ 2}(\overline{\Omega}\times {\mathbb{R}})$.
For all $(x,t)\in\Omega\times {\mathbb{R}}$,  
$0 <\alpha_i\leq a_i(x,t)\leq A_i$,  
$0<\beta_i\leq b_i(x,t)\leq B_i$,
$ 0<\gamma_i\leq c_i(x,t)\leq C_i$, where
$\alpha_i, \beta_i, \gamma_i,  A_i, B_i, C_i \in {\mathbb{R}}$.
Observe that under competition, the growth of each population is
reduced at a rate proportional to the size of the population of its
competitor. We are concerned with the existence of positive
solutions for the problem \eqref{Syexam3}, which can be translated
as the co-existence
of the two populations without  asymptotic extinction.

Note that the reaction functions are quasimonotone nonincreasing in  
$[0,\infty)\times[0,\infty)$.
To apply Theorem \ref{Syt1}, we need to have an ordered
sub-solution $(\underline{u},\underline{v})$ and  super-solution
$(\overline{u},\overline{v})$ of \eqref{Syexam3}. As in Definition
\ref{defsubsol}, they need to satisfy the following inequalities
\begin{equation}
\begin{gathered} \label{compet1}
\frac{\partial \underline{u}}{\partial t}-\triangle
\underline{u} \le \underline{u}(a_1(x,t)-b_1(x,t)
\underline{u}-c_1(x,t) \overline{v}) \quad\text{in } \Omega\times
{\mathbb{R}},\\
\frac{\partial \overline{v}}{\partial t}-\triangle \overline{v}
\geq \overline{v}(a_2(x,t)-b_2(x,t) \overline{v}-c_2(x,t)
\underline{u})\quad\text{in } \Omega\times
{\mathbb{R}},\\
 \frac{\partial \overline{u}}{\partial t}-\triangle
\overline{u} \ge \overline{u}(a_1(x,t)-b_1(x,t)
\overline{u}-c_1(x,t) \underline{v}) \quad\text{in } \Omega\times
{\mathbb{R}},\\
\frac{\partial \underline{v}}{\partial t}-\triangle \underline{v}
\leq \underline{v}(a_2(x,t)-b_2(x,t) \underline{v}-c_2(x,t)
\overline{u})\quad\text{in } \Omega\times
{\mathbb{R}},\\
\frac{\partial \underline{u}}{\partial\nu}
 \le 0\le\frac{\partial \overline{u}}{\partial\nu} 
 \quad\text{on } \partial{\Omega}\times {\mathbb{R}},\\
\frac{\partial \underline{v}}{\partial\nu} 
 \leq 0\leq \frac{\partial \overline{v}}{\partial\nu}\quad\text{on } 
\partial{\Omega}\times {\mathbb{R}},\\
 \sup_{\Omega\times
{\mathbb{R}}}\{|\underline{u}|,|\overline{u}|,|\underline{v}|,
|\overline{v}|\}<\infty.
\end{gathered}
\end{equation}
Choose $(\overline{u},\overline{v})=(M_1, M_2)$ such that 
$M_i\geq \frac{A_i}{\beta_i}\; (i=1,2)$. Pick
$(\underline{u},\underline{v})=(\varepsilon_1,\varepsilon_2)$ such
that $\varepsilon_i\leq \frac{\alpha_i\beta_j-A_jC_i}{\beta_jB_i}$
with $C_i<\frac{\alpha_i\beta_j}{A_j}$ for $i,j=1,2$ and 
($i\neq j)$.  With this choice  of sub- and super-solutions,  one can see
that $(\underline{u},\underline{v})\leq \overline{u},\overline{v})$
and inequalities  in \eqref{compet1}  are satisfied. Therefore, it
follows from Theorem \ref{Syt1} that problem  \ref{Syexam3} has a
positive solution $(u,v)$ such that 
$(\varepsilon_1, \varepsilon_2 )\leq (u,v)\leq (M_1, M_2)$
in $\overline{\Omega}\times {\mathbb{R}}$. Thus, $(u,v)$ does not tend
to zero as $t$ tends to $\pm \infty$, for each $x\in\Omega$.
\end{example}

\begin{example} \rm (Competitive Model with Nonlinear Boundary Conditions)
\begin{equation}\label{Syexam4}
\begin{gathered}
 \frac{\partial u}{\partial t}-\triangle u =
u(a_1(x,t)-b_1(x,t) u-c_1(x,t) v) \quad\text{in } \Omega\times
{\mathbb{R}},\\
\frac{\partial v}{\partial t}-\triangle v = v(a_2(x,t)-b_2(x,t)
v-c_2(x,t) u)\quad\text{in } \Omega\times
{\mathbb{R}},\\
\frac{\partial u}{\partial\nu}= u^m(\delta_1-u-\sigma_1v)\quad\text{on } 
 \partial{\Omega}\times {\mathbb{R}},\\
\frac{\partial v}{\partial\nu}=v^m(\delta_2-v-\sigma_2u) \quad\text{on }
  \partial{\Omega}\times {\mathbb{R}},\\
 \sup_{\Omega\times
{\mathbb{R}}}\{|u|,|v|\}<\infty.
\end{gathered}
\end{equation}
where $0<\delta_i, \sigma_i\in {\mathbb{R}}$,
$a_i:\Omega\times {\mathbb{R}}\to {\mathbb{R}}$,
$b_i:\Omega\times {\mathbb{R}}\to {\mathbb{R}}$,
$c_i:\Omega\times {\mathbb{R}}\to {\mathbb{R}}$ are in
$C_{\rm{loc}}^{\mu,\mu/ 2}(\overline{\Omega}\times {\mathbb{R}})$.
For all $(x,t)\in\Omega\times {\mathbb{R}}$, 
$ 0 <\alpha_i\leq a_i(x,t)\leq A_i,  0<\beta_i\leq b_i(x,t)\leq B_i$,
$ 0<\gamma_i\leq c_i(x,t)\leq C_i$, where
$\alpha_i, \beta_i, \gamma_i,  A_i, B_i, C_i \in {\mathbb{R}}$.
\end{example}

Note that the nonlinearities satisfy the quasimonotone nondecreasing
property. Choose $(\overline{u},\overline{v})=(M_1,M_2)$ with
$M_i\ge \delta_i+\frac{A_i}{\beta_i} (i=1,2)$, and pick
$(\underline{u},\underline{v})=(\varepsilon_1,\varepsilon_2)$, where
$0<\varepsilon_i<\min\{\delta_i-\sigma_iM_j, \;\alpha_i-C_iM_j\}$
 with
$C_i<\frac{\alpha_j}{\delta_j+\frac{A_j}{\beta_j}}$ 
($i,j=1,2$ with $i\neq j$). Then by Theorem \ref{Syt1}, the system
\eqref{Syexam4} has a positive solution $(u,v)$ such that
$(\varepsilon_1, \varepsilon_2)\leq (u,v)\leq (M_1, M_2)$
in $\overline{\Omega}\times {\mathbb{R}}$. Thus, $(u,v)$ does
not tend to zero as $t$ tends to $\pm \infty$.


\begin{example}\rm 
(Nonlinearities with no one-sided Lipschitz conditions)
\begin{equation}\label{Syexam5}
\begin{gathered}
 \frac{\partial u}{\partial t}-\triangle u =
f_1(x,t,u,v)\quad\text{in } \Omega\times {\mathbb{R}},\\
\frac{\partial v}{\partial t}-\triangle v =f_2(x,t,u,v)\quad\text{in }
 \Omega\times {\mathbb{R}},\\
\frac{\partial u}{\partial\nu}= 0=\frac{\partial v}{\partial\nu}\quad\text{on } 
\partial{\Omega}\times {\mathbb{R}},\\
 \sup_{\Omega\times {\mathbb{R}}}\{|u|,|v|\}<\infty.
\end{gathered}
\end{equation}
where  \begin{equation*}
f_1(x,t,u,v)=\begin{cases} -b_1(x,t)u^{\mu} v& \; \text{if } 0\le u, v<\infty, \,{\rm for \;some } \, 0<\mu<1\\
0 &  \; \text{if } -\infty\le u, v\leq 0
\end{cases}
\end{equation*}
where  
\begin{equation*}
f_2(x,t,u,v)=\begin{cases} 
-b_2(x,t)v^{\mu} u &  \text{if } 0\le u, v<\infty, \text{ for some } 0<\mu<1\\
0 &  \text{if } -\infty\le u, v\leq 0
\end{cases}
\end{equation*}
and $b_i\in L^{\infty}({\Omega\times {\mathbb{R}}})$
such that $0<\beta\le b_i(x,t)\le B$ for a.e 
$(x,t)\in \Omega\times {\mathbb{R}}$, where $\beta, B\in \mathbb{R}$.


 It is easily seen that $(\underline{u},\underline{v})=(0,0)$ and 
$(\overline{u},\overline{v})=(K,K)$ (where $0<K\in \mathbb{R}$) 
are ordered subsolution and supersolution of problem \eqref{Syexam5}. 
Therefore, by Theorem \ref{Syt1}, there exists a solution $(u,v)$ 
of problem \eqref{Syexam5} such that $(0,0)\leq (u,v)\leq (K,K)$.
\end{example}

\begin{thebibliography}{9}
\bibitem{HA71} H. Amann;
\emph{On the existence of positve solutions of nonlinear
elliptic boundary value problems,} Indiana Univ. Math. J. 21 (1971),
125-146.

\bibitem{HA76S} H. Amann;
 \emph{Nonlinear Elliptic Equations with Nonlinear Boundary
Conditions,} Proc. 2nd Schrveningen Conf. Diff. Eq., North-Holland
Math. Studies 21, pp. 43-64, 1976.

\bibitem{HA88} H. Amann;
\emph{Dynamic theory of quasilinear parabolic equations}, 
I. Abstract evolution equations, Nonlinear Anal. 12 (1988), 895–919.

\bibitem{CC03S} R.  S. Cantrell and  C. Cosner, \emph{Spatial Ecology Via Reaction-Diffusion
Equations,} John Wiley\& Sons Ltd, England, 2003.

\bibitem{CC09} Z. Chen, Y. Cui;
 \emph{Existence and Uniqueness of solutions to a semilinear elliptic system},
 Electronic Journal of Differential Equations, Vol. 2009(2009), No. 134, pp. 1–7.

\bibitem{DO00} C. De Coster, P. Omari;
\emph{Unstable periodic solutions of a parabolic problem in the presence 
of non-well-ordered lower and upper solutions,} J. Func. Anal. 175 (2000), 52-88.

\bibitem{PF64S} P. Fife;
 \emph{Solutions of Parabolic Boundary Problems Existing for All
Time,} Arch. Rational Mech. Anal.  16 (1964),  155 - 186.

\bibitem{PF79S} P. Fife;
\emph{Mathematical Aspects of Reacting and Diffusing
Systems,} Lecture Notes in Biomathematics 28, Springer-Verlag, New
York, 1979.

\bibitem{FM90} D. G. de Figueiredo, E. Mitidieri;
 \emph{Maximum Principles for Cooperative Elliptic Systems},
 Comptes Rendus Acad.Sc.Paris, 310, (1990), 49-52.

\bibitem{AF64S} A. Freidman;
\emph{Partial Differenial Equations of Parabolic Type},
Prentice-hall, NJ, 1964.

\bibitem{PH91S} P. Hess;
\emph{Periodic-parabolic Boundary Value Problems and
Positivity,} Longman Group, UK, 1991.

\bibitem{NVK96S} N. V. Krylov;
\emph{Lectures on Elliptic and Parabolic Equations in
H\"{o}lder Spaces}, Graduate Studies in Mathematics 12, AMS,
Providence, Rhodes Island, 1996.

\bibitem{LSU68S} O. A. Lady\v{z}enskaja, V. A. Solonnikov,  N. N.
Ural'ceva;
\emph{Linear and Quasilinear Equations of Parabolic
Type}, AMS translations of Mathematical Monographs, Vol. 23,
Providence, Rhodes Island, 1968.

\bibitem{AL82S} A. Leung;
\emph{A Semilinear Reaction-Diffusion Prey-Predator System
with Nonlinear Coupled Boundary Conditions: Equilibrium and
Stability}, Indiana University Mathematics Journal, Vol. 31, No. 2,
1982.

\bibitem{NM11} N. Mavinga and M. N. Nkashama;
 \emph{Strong full bounded solutions of nonlinear Parabolic Equations 
with nonlinear boundary conditions,}  
Nonlinear Analysis {{74}} (2011),  pp. 5171--5188.

\bibitem{NM10}  N. Mavinga, M. N. Nkashama;
 \emph{Bounded Solutions of Nonlinear Parabolic Equations,}  
Electronic Journal of Differential Equations, Conference {{19}} (2010),  
pp. 197--205.

\bibitem{WM60} W. Mlak;
\emph{Parabolic differential inequalities and
Chaplighin's method,} Ann. Polon. Math. 8 (1960), 139-152.

\bibitem{MP09} T. Moussai, R. Precup;
\emph{Existence of solutions for second-order differential equations and 
systems on inifinite intervals}, 
Electronic Journal of Differential Equations, Vol. 2009(2009), No. 94, pp. 1–13.

\bibitem{P92S} C. V. Pao;
\emph{Nonlinear Parabolic and Elliptic Equations,} Plenum Press, New York, 1992.

\bibitem{PW67S} M. H. Protter, H. F. Weinberger;
\emph{Maximum Principles in Differential equations,} 
Prentice-Hall, Englewood Cliffs, NJ, 1967.

\bibitem{DH73S} D. H. Sattinger;
\emph{Topics in Stability and Bifurcation Theory,}  
Lecture Notes in Mathematics, No. 309, Springer-Verlag,
New York, 1979.

\bibitem{GS71S} G. R. Sell;
\emph{Topological dynamics and ordinary differential equations,} 
Van Nostrand Reinhold Compagny, 1971.

\bibitem{SZ94} J. Serrin, H. Zou;
\emph{ Non-existence of positive solutions of semilinear elliptic systems}, 
Discourses in Mathematics and its Applications, no. 3, Department of Mathematics,
 Texas A\&M University, College Station, Texas, (1994), 55-68.

\bibitem{SY98S} W. Shen, Y. Yi;
\emph{Convergence in almost periodic Fisher and
Kolmogorov models,} J. Math. Biol. 37,  84-102, 1998.

\bibitem{HS95} H. Smith;
\emph{Monotone Dynamical Systems,} Mathematical Surveys and Monographs 41, 
AMS, Providence, RI, 1995.

\end{thebibliography}

\end{document}