We study the following system of nonlinear Schrödinger equations: $$\begin{aligned} \left\{ \begin{array}{l} -\varepsilon ^2\Delta u +a(x) u = f(u)+\lambda v, \quad x\in \mathbb R ^N, \\ -\varepsilon ^2\Delta v +b(x) v =g(v)+\lambda u, \quad x\in \mathbb R ^N,\\ u,v >0 \,\,\,\hbox {in}\,\,\,\mathbb R ^N,\quad u, v \in H^1 (\mathbb R ^N), \end{array}\right. \end{aligned}$$ - ε 2 Δ u + a ( x ) u = f ( u ) + λ v , x ∈ R N , - ε 2 Δ v + b ( x ) v = g ( v ) + λ u , x ∈ R N , u , v > 0 in R N , u , v ∈ H 1 ( R N ) , where $$N\ge 3$$ N ≥ 3 , $$\varepsilon , \lambda >0$$ ε , λ > 0 , and $$a, b, f, g$$ a , b , f , g are continuous functions. Under very general assumptions on both the potentials $$a, b$$ a , b and the nonlinearities $$f, g$$ f , g , for small $$\lambda >0$$ λ > 0 and $$\varepsilon >0$$ ε > 0 , we obtain positive solutions of this coupled system via pure variational methods. The asymptotic behaviors of these solutions are also studied either as $$\varepsilon \rightarrow 0$$ ε → 0 or as $$\lambda \rightarrow 0$$ λ → 0 .