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