We consider the following convective Neumann systems: \( ( \mathrm{S} ) \quad \textstyle\begin{cases} -\Delta _{p_{1}}u_{1}+ \frac{ \vert \nabla u_{1} \vert ^{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{ \vert \nabla u_{2} \vert ^{p_{2}}}{u_{2}+\delta _{2} }=f_{2}(x,u_{1},u_{2}, \nabla u_{1},\nabla u_{2}) & \text{in } \Omega , \\ \vert \nabla u_{1} \vert ^{p_{1}-2}\frac{\partial u_{1}}{\partial \eta }=0= \vert \nabla u_{2} \vert ^{p_{2}-2}\frac{\partial u_{2}}{\partial \eta } & \text{on } \partial \Omega ,\end{cases} \) where Ω is a bounded domain in $\mathbb{R}^{N}$ ( $N\geq 2$ ) with a smooth boundary ∂Ω, $\delta _{1}, \delta _{2} >0$ are small parameters, η is the outward unit vector normal to ∂Ω, $f_{1}, f_{2}:\Omega \times \mathbb{R}^{2}\times \mathbb{R}^{2N} \rightarrow \mathbb{R}$ are Carathéodory functions that satisfy certain growth conditions, and $\Delta _{p_{i}}$ ( $1< p_{i}< N$ , $i=1,2$ ) are the p-Laplace operators $\Delta _{p_{i}}u_{i}=\operatorname{div}(|\nabla u_{i}|^{p_{i}-2}\nabla u_{i})$ for $u_{i}\in W^{1,p_{i}}(\Omega )$ . To prove the existence of solutions to such systems, we use a subsupersolution method. We also obtain nodal solutions by constructing appropriate subsolution and supersolution pairs. To the best of our knowledge, such systems have not been studied yet.
