Multiple solutions for a Neumann system involving subquadratic nonlinearities

In this paper we consider the model semilinear Neumann system $$\{\begin{array}{lll}-\mathrm{\Delta }u+a(x)u=\lambda c(x)F_u(u,v)& \mathrm{i}\mathrm{n}& \mathrm{\Omega },\\ -\mathrm{\Delta }v+b(x)v=\lambda c(x)F_v(u,v)& \mathrm{i}\mathrm{n}& \mathrm{\Omega },\\ \frac{\mathrm{\partial }u}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }v}{\mathrm{\partial }\nu }=0& \mathrm{o}\mathrm{n}& \mathrm{\partial }\mathrm{\Omega },\end{array}$$ where $\mathrm{\Omega }\subset {\mathbb{R}}^N$ is a smooth open bounded domain, $\nu$ denotes the outward unit normal to $\mathrm{\partial }\mathrm{\Omega }$, $\lambda \ge 0$ is a parameter, $a,b,c\in L_{+}^{\mathrm{\infty }}(\mathrm{\Omega })\setminus \{0\},$ and $F\in C^1({\mathbb{R}}^2,\mathbb{R})\setminus \{0\}$ is a nonnegative function which is subquadratic at infinity. Two nearby numbers are determined in explicit forms, $\underset{‾}{\lambda }$ and $\bar{\lambda }$ with $0<\underset{‾}{\lambda }\le \bar{\lambda }$, such that for every $0\leq \lambda<\underline \lambda$,system$(N_\lambda)$ has only the trivial pair of solution, while for every $\lambda >\bar{\lambda }$, system $(N_{\lambda })$ has at least two distinct nonzero pairs of solutions.