  We consider the following convective Neumann
systems:\begin{equation*}\left(\mathrm{S}\right)\qquad\left\{\begin{array}{ll}-\Delta_{p_1}u_1+\frac{|\nabla
u_1|^{p_1}}{u_1+\delta_1}=f_1(x,u_1,u_2,\nabla u_1,\nabla u_2) &
\text{in}\;\Omega,\\ -\Delta _{p_2}u_2+\frac{|\nabla
u_2|^{p_2}}{u_2+\delta_2}=f_2(x,u_1,u_2,\nabla u_1,\nabla
u_2)&\text{in}\;\Omega, \\ |\nabla u_1|^{p_1-2}\frac{\partial u_1}{\partial
\eta }=0=|\nabla u_2|^{p_2-2}\frac{\partial u_2}{\partial
\eta}&\text{on}\;\partial\Omega,\end{array}\right.\end{equation*}where $\Omega$
is a bounded domain in $\mathbb{R}^{N}$ ($N\geq 2$) with a smooth boundary
$\partial\Omega$,$\delta_1,\,\delta_2 >0$ are small parameters, $\eta$ is the
outward unit vector normal to $\partial \Omega,$
$f_1,\,f_2:\Omega\times\mathbb{R}^2\times\mathbb{R}^{2N}\rightarrow \mathbb{R}$
are Carath\'{e}odory functions that satisfy certain growth conditions, and
$\Delta _{p_i}$ ($1<p_i<N,$ for $i=1,2$) are the $p$-Laplace operators $\Delta
_{p_i}u_i=\mathrm{div}(|\nabla u_i|^{p_i-2}\nabla u_i)$,for every $\,u_i\in
W^{1,p_i}(\Omega).$ In order to prove the existence of solutions to such
systems, we use a sub-supersolution method. We also obtain nodal solutions by
constructing appropriate sub-solution and super-solution pairs. To the best of
our knowledge, such systems have not been studied yet.
