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.