We consider the problem \begin{equation*} \begin{cases} \varepsilon^2 \Delta u - u + f(u) = 0, u > 0 & \mbox{in } \Omega\\ \frac{\partial u}{\partial \nu} = 0 & \mbox{on } \partial\Omega, \end{cases} \end{equation*} where $\Omega$ is a bounded smooth domain in $R^N$, $\ve>0$ is a small parameter and $f$ is a superlinear, subcritical nonlinearity. It is known that this equation possesses multiple boundary spike solutions that concentrate, as $\epsilon$ approaches zero, at multiple critical points of the mean curvature function $H(P)$, $P \in \partial \Omega$. It is also proved that this equation has multiple interior spike solutions which concentrate, as $\ep\to 0$, at {\it sphere packing\/} points in $\Om$. In this paper, we prove the existence of solutions with multiple spikes {\it both\/} on the boundary and in the interior. The main difficulty lies in the fact that the boundary spikes and the interior spikes usually have different scales of error estimation. We have to choose a special set of boundary spikes to match the scale of the interior spikes in a variational approach.