In this paper, we consider nonlinear Schrödinger equations of the following type:−Δu(x)+ V(x)u(x) − q(x)|u(x)|σu(x) = λu(x), x ∈ ℝN, u ∈ H1(ℝN)\{0},where N ≥ 2 and σ > 0. We concentrate on situations where the potential function V appearing in the linear part of the equation is of Coulomb type; by this we mean potentials where the spectrum of the linear operator −Δ + V consists of an increasing sequence of eigenvalues λ1, λ2,… followed by an interval belonging to the essential spectrum.We study, for λ kept fixed inside a spectral gap or below λ1, the existence of multiple solution pairs, as well as the bifurcation behaviour of these solutions when λ approaches a point of the spectrum from the left-hand side. Our method proceeds by an analysis of critical points of the corresponding energy functional. To this end, we derive a new variational characterization of critical levelsc0 (λ) ≤ c1(λ) ≤ c2(λ) ≤ ⋯ corresponding to an infinite set of critical points.We derive such a multiplicity result even if some of the critical values cn(λ) coincide; this seems to be a major advantage of our approach. Moreover, the characterization of these values cn(λ) is suitable for an analysis of the bifurcation behaviour of the corresponding generalized solutions.The approach presented here is generic; for instance, it can be applied when V and q are periodic functions. Such generalizations are briefly described in this paper and will be the object of a forthcoming article.