The existence and multiplicity of weak solutions is established for a class of concave–convex elliptic systems of the form: { − Δ p u + m ( x ) | u | p − 2 u = λ a ( x ) | u | γ − 2 u , x ∈ Ω , − Δ p v + m ( x ) | v | p − 2 v = μ b ( x ) | v | γ − 2 v , x ∈ Ω , | ∇ u | p − 2 ∂ u ∂ n = α α + β | u | α − 2 u | v | β , | ∇ v | p − 2 ∂ v ∂ n = β α + β | u | α | v | β − 2 v , x ∈ ∂ Ω . Here Δ p denotes the p -Laplacian operator defined by Δ p z = div ( | ∇ z | p − 2 ∇ z ) , p > 2 , Ω ⊂ R N is a bounded domain with smooth boundary, α > 1 , β > 1 , 2 < α + β < p < γ < p ∗ ( p ∗ = p N N − p if N > p , p ∗ = ∞ if N ≤ p ), ∂ ∂ n is the outer normal derivative, ( λ , μ ) ∈ R 2 ∖ { ( 0 , 0 ) } , the weight m ( x ) is a positive bounded function and a ( x ) , b ( x ) ∈ C ( Ω ) are functions which change sign in Ω . Our technical approach is based on the Nehari manifold which is similar to the fibering method of Drabek and Pohozaev (1997) [29] together with the recent idea from Brown and Wu (2008) [10].