We consider the problem { ε 2 Δ u − u + f ( u ) = 0 , u > 0 , in Ω , ∂ u ∂ v = 0 on ∂ , where Ω is a bounded smooth domain in R N , ɛ > 0 is a small parameter, and ƒ is a superlinear, subcritical nonlinearity. It is known that this equation possesses boundary spike solutions that concentrate, as ɛ approaches zero, at a critical point of the mean curvature function H(P), P ε ∂Ω. It is also proved that this equation has single interior spike solutions at a local maximum point of the distance function d(P, ∂Ω), P ε Ω . In this paper, we prove the existence of interior K −peak (K ⩾ 2) solutions at the local maximum points of the following function φ(P 1 , P 2 ,…, p k ) = min i, k, l = 1, …, K; k ≠ 1 (d(P i , ∂Ω), 1/2 ¦P k −p l ¦) . We first use the Liapunov-Schmidt reduction method to reduce the problem to a finite dimensional problem. Then we use a maximizing procedure to obtain multiple interior spikes. The function φ(P 1 , …, p k ) appears naturally in the asymptotic expansion of the energy functional.