  We consider the following elliptic system with Neumann boundary:
\begin{equation} \begin{cases} -\Delta u + \mu u=v^p, &\hbox{in } \Omega,
\\-\Delta v + \mu v=u^q, &\hbox{in } \Omega, \\\frac{\partial u}{\partial n} =
\frac{\partial v}{\partial n} = 0, &\hbox{on } \partial\Omega, \\u>0,v>0,
&\hbox{in } \Omega, \end{cases} \end{equation} where $\Omega \subset
\mathbb{R}^N$ is a smooth bounded domain, $\mu$ is a positive constant and
$(p,q)$ lies in the critical hyperbola: $$ \dfrac{1}{p+1} + \dfrac{1}{q+1}
=\dfrac{N-2}{N}. $$ By using the Lyapunov-Schmidt reduction technique, we
establish the existence of infinitely many solutions to above system. These
solutions have multiple peaks that are located on the boundary $\partial
\Omega$. Our results show that the geometry of the boundary $\partial\Omega,$
especially its mean curvature, plays a crucial role on the existence and the
behaviour of the solutions to the problem.
