We consider the following singularly perturbed Neumann problem: epsilon(2) Delta u - u + f(u) = 0 in Omega, u > 0 in Omega, (partial derivative u)/(partial derivative v) = 0 on partial derivative Omega, where Delta = Sigma(N)(i=1) partial derivative(2)/partial derivative x(i)(2) is the Laplace operator, epsilon > 0 is a constant, Omega is a bounded, smooth domain in RN with its unit outward normal v, and f is superlinear and subcritical. A typical f is f(u) = u(p) where 1 < p < +infinity when N = 2 and 1 < p < (N+2)/(N-2)when N >= 3. We show that there exists an epsilon(0) > 0 such that for 0 < epsilon < epsilon(0) and for each integer K bounded by 1 <= K <= alpha(N,Omega,f)/epsilon(N)(vertical bar ln epsilon vertical bar)(N) where alpha(N,Omega,f) is a constant depending on N, Omega, and f only, there exists a solution with K interior peaks. (An explicit formula for alpha(N,Omega,f) is also given.) As a consequence, we obtain that for c sufficiently small, there exists at least [alpha(N,Omega f)/epsilon(N)(vertical bar ln epsilon vertical bar)(N)] number of solutions. Moreover, for each m is an element of (0, N) there exist solutions with energies in the order of epsilon(N-m). (c) 2006 Wiley Periodicals, Inc.