  In this article we shall study the following elliptic system with
coefficients:
  \begin{equation}\notag
  \left\{\begin{aligned}
  -\epsilon^2\Delta u +c(x)u=b(x)|v|^{q-1}v, &\text{ and } -\epsilon^2\Delta v
+c(x)v=a(x) |u|^{p-1}u &&\text{in } \Omega \newline
  u>0, \ v>0 \text{ in } \Omega, &\text{ and }\quad\frac{\partial
u}{\partial\nu} = 0 = \frac{\partial v}{\partial\nu} &&\text{on }\partial\Omega
  \end{aligned}
  \right.
  \end{equation} where $\Omega$ is a smooth bounded domain in $R^n, n\geq 3$.
The coefficients $a(x), b(x)$ and $c(x)$ are positive bounded smooth functions.
We shall study the existence of point concentrating solutions and discuss the
role of the coefficients to determine the concentration profile of the
solutions. We have also discussed some applications of our main theorem towards
the existence of solutions concentrating on higher-dimensional orbits.
