In this paper, we study the following linearly coupled elliptic system (Pε) \(\begin{cases}-\varepsilon^{2}\Delta u+P(x)u=u^{3}+\lambda(x)v & \text{in} \ \Omega,\\ -\varepsilon^{2}\Delta v+Q(x)v=v^{3}+\lambda(x)u & \text{in} \ \Omega,\\ u>0,\ v>0 & \text{in} \ \Omega,\\ {\partial u\over \partial n}={\partial v\over \partial n}=0 & \text{on}\ \partial \Omega,\end{cases}\) where ε > 0, Ω is smooth and bounded in ℝ3 with boundary ∂Ω, and n is the outer normal vector defined on ∂Ω.
Let ω be the unique positive radial solution of the well-known equation \(-\Delta \omega+\omega=\omega^{3}, \quad\omega \in H^{1}(\mathbb{R}^{3}) ,\) and μ1 < 0 be the first eigenvalue of the operator −Δ + id − 3w2 defined on H1(ℝ3). Assume that \(P(x), \ Q(x), \ \lambda(x) \in C^{1}(\overline{\Omega})\) satisfy 0 < λ(x) < min{P(x), Q(x)}, \(P(x)=Q(x)=a_{i}> 0, \quad \lambda(x)=\lambda_{i}\in (0,a_{i}),\quad \forall x \in N_{i}, \ i=1,2,\ldots, K,\) where (ai − λi)μ1 + 2λi ≠ 0, {Ni ⊂ ∂Ω∣i = 1, 2,…, K} are pair-wise disjoint neighborhoods of the local minima (maxima) of the mean curvature H(P), P ∈ ∂Ω. Via Lyapunov–Schmidt reduction method, we may construct a solution with K peaks to the system (Pε) with each peak being on ∂Ω and locating near these local minima (maxima) points.
