Multiple solutions for a Neumann system involving subquadratic nonlinearities

In this paper we consider the model semilinear Neumann system $$\left\{ \begin{array}{lll} -\Delta u+a(x)u=\lambda c(x) F_u(u,v)& {\rm in} & \Omega,\\ -\Delta v+b(x)v=\lambda c(x) F_v(u,v)& {\rm in} & \Omega,\\ \frac{\partial u}{\partial \nu}=\frac{\partial v}{\partial \nu}=0 & {\rm on} & \partial\Omega, \end{array}\right.$$ where $\Omega\subset \mathbb R^N$ is a smooth open bounded domain, $\nu$ denotes the outward unit normal to $\partial \Omega$, $\lambda\geq 0$ is a parameter, $a,b,c\in L_+^\infty(\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, $\underline \lambda$ and $\overline \lambda$ with $0<\underline\lambda\leq \overline \lambda$, such that for every $0\leq \lambda<\underline \lambda$, system $(N_\lambda)$ has only the trivial pair of solution, while for every $\lambda>\overline \lambda$, system $(N_\lambda)$ has at least two distinct nonzero pairs of solutions.