Nodal solutions for Neumann systems with gradient dependence

We consider the following convective Neumann systems: ( S ) { − Δ p 1 u 1 + | ∇ u 1 | p 1 u 1 + δ 1 = f 1 ( x , u 1 , u 2 , ∇ u 1 , ∇ u 2 ) in Ω , − Δ p 2 u 2 + | ∇ u 2 | p 2 u 2 + δ 2 = f 2 ( x , u 1 , u 2 , ∇ u 1 , ∇ u 2 ) in Ω , | ∇ u 1 | p 1 − 2 ∂ u 1 ∂ η = 0 = | ∇ u 2 | p 2 − 2 ∂ u 2 ∂ η on ∂ Ω , where Ω is a bounded domain in R N ( N ≥ 2 ) with a smooth boundary ∂ Ω, δ 1 , δ 2 > 0 are small parameters, η is the outward unit vector normal to ∂ Ω, f 1 , f 2 : Ω × R 2 × R 2 N → R are Carathéodory functions that satisfy certain growth conditions, and Δ p i ( 1 < p i < N , i = 1 , 2 ) are the p -Laplace operators Δ p i u i = div ( | ∇ u i | p i − 2 ∇ u i ) for u i ∈ W 1 , p i ( Ω ) . 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.