In this paper, we consider nonlinear Schrodinger equations of the following type: -Delta u(x) + V (x) u(x) -q(x)|u(x)|sigma u(x) =lambda u(x), x is an element of R-N, u is an element of H-1(R-N) \ {0}, where N >= 2 and sigma > 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 -Delta + V consists of an increasing sequence of eigenvalues lambda(1), lambda(2), followed by an interval belonging to the essential spectrum. We study, for lambda kept fixed inside a spectral gap or below lambda(1), the existence of multiple solution pairs, as well as the bifurcation behaviour of these solutions when lambda 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 levels c(0)(lambda) <= c(1)(lambda) <= c(2)(lambda) <= center dot center dot center dot corresponding to an infinite set of critical points. We derive such a multiplicity result even if some of the critical values c(n)(lambda) coincide; this seems to be a major advantage of our approach. Moreover, the characterization of these values c(n)(lambda) 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.