Motivated by problems arising in nonlinear optics and Bose-Einstein condensates, we consider in R-N, with N <= 3, the following system of coupled Schrodinger equations: {-Delta u(i) + lambda V-i(x)u(i) = u(i)Sigma(d)(j=1)alpha(ij)u(j)(2), i=1, ..., d, u(i)>= 0, lim(vertical bar x vertical bar ->infinity) u(i)(x) = 0, where lambda > 0 is a parameter. alpha(ij) = alpha(ji) are positive constants, and V-i non-negative given potentials. We assume that the interior of boolean AND(d)(i=1) V-i(-1) (0) admits ni connected components Omega(1), ..., Omega(m) which are of class C-1, and isolated in each V-i(-1) (0). For each non-empty J subset of {1, ..., m}, we prove that the system admits for any lambda large a multi-bump solution u(lambda) :R-N -> R-d which is small in R-N\boolean OR(j is an element of J)Omega(j), and on each Omega(j) (j is an element of J) close in H-1-norm to a least energy solution of the limit problem: -Delta u(i) = u(i)Sigma(d)(j=1)alpha(ij)u(j)(2), i=1, ..., d, subjected to homogeneous Dirichlet boundary condition. An explicit condition on the matrix (alpha(ij)) is given to ensure our solutions have at least two positive components. (C) 2011 Elsevier Inc. All rights reserved.