Summary: We consider the following convective Neumann systems: \[ (\text{S}) \quad \ \begin{cases} -\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{cases} \] 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} \to \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.