The goal of this paper is to study the existence of non-trivial weak solution for the following system of nonlinear elliptic equations: \begin{eqnarray*} &&-\mathrm{div}(h_1(x)\nabla u)+a(x)u=f(x,u,v) \text{ in } \Omega\\ && -\mathrm{div}(h_2(x)\nabla v)+b(x)v=g(x,u,v) \text{ in } \Omega \end{eqnarray*} with Neumann condition: \begin{gather*} \dfrac{\partial u}{\partial n}=0,\quad \dfrac{\partial v}{\partial n}=0\\ u(x)\longrightarrow 0, \quad v(x)\longrightarrow 0 \text{ as } |x|\longrightarrow +\infty \end{gather*} where $\Omega \subset R^N, N \ge 3$ is an unbounded domain with smooth bounded boundary $\partial\Omega$, and $h_i(x)\in L^1_{loc}(\overline{\Omega})$, $i=1,2$, $\overline{\Omega}=\Omega\cup\partial\Omega$. The solutions will be obtained in a subspace of the space $H^1(\Omega)$ and the proofs rely essentially on a variation of the Mountain Pass Theorem in [7].