  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.
