In this paper, we study the following nonlinear Hartree system \(\left\{ \begin{array}{ll} -\Delta u +P(x)u=\mu _1 \phi _u u+\beta \phi _vu,~~& x\in {\mathbb {R}}^3, \\ -\Delta v +Q(x)v=\mu _2 \phi _v v+\beta \phi _uv,~~& x\in {\mathbb {R}}^3,\\ \end{array} \right. \) where \( \phi _u(x):=\int _{{\mathbb {R}}^3}\frac{u^2(y)}{|x-y|}dy\) for any \(u\in H^1({\mathbb {R}}^3),\) and P(x) and Q(x) are continuous bounded radial functions. When P(x) and Q(x) satisfy some properties of decay at infinity and \(\beta \) takes some appropriate value, we construct infinitely many both synchronized and segregated solutions by the finite dimensional reduction argument. Moreover, when k is a sufficiently large even integer, we can also construct the positive solutions with k peaks and sign-changing solutions with k peaks under the same conditions. During the process of constructing the synchronized solutions, our results contain the cases \(\mu _1\le 0\) or \(\mu _2\le 0\) , which correspond to the Poisson problem and are new in this aspect. For the segregated solutions, it is very surprising that different decay rates of the potentials cause different ranges for \(\beta \) . Our results can be seen as extending that of [29](Peng and Wang, ARMA, 2013, 208:305-339) to the non-local problem and those of [16](Hu, Jevnikar and Xie, Commun. Contemp. Math., 2023, 25:Paper No. 2350008, 19 pp) to the system. However, studying the above problem, we encounter some new difficulties, caused by the non-local terms, since there is few knowledge about the non-degeneracy of the positive solutions to the limit system and the error estimates of the coupled non-local term in constructing segregated solutions. Also, we obtain some results that differ from those in the local case.
