We study the following system of nonlinear Schrodinger equations: { -epsilon(2)Delta u + a(x)u = f(u) + lambda v, x is an element of R-N, -epsilon(2)Delta v + b(x)v = g(v) + lambda u, x is an element of R-N, u,v > 0 in R-N , u, v is an element of H-1 (R-N), where N >= 3, epsilon, lambda > 0, and a, b, f, g are continuous functions. Under very general assumptions on both the potentials a, b and the nonlinearities f, g, for small lambda > 0 and epsilon > 0, we obtain positive solutions of this coupled system via pure variational methods. The asymptotic behaviors of these solutions are also studied either as epsilon -> 0 or as lambda -> 0 .