Bound states of the generalized spiked harmonic oscillator potential are calculated accurately by using the generalized pseudospectral method. Energy eigenvalues, various expectation values, radial densities are obtained through a nonuniform, optimal spatial discretization of the radial Schrödinger equation efficiently. Ground and excited states corresponding to arbitrary values of n and ℓ are reported for potential parameters covering a wide range of interaction. Both weak and strong coupling is considered. The effect of potential parameters on eigenvalues and densities are discussed. Almost all of the results are reported here for the first time, which could be useful for future studies.The spiked harmonic oscillator (SHO) defined by the Hamiltonian H = −d 2 /dr 2 + r 2 + λr −α , where r ∈ [0, ∞] and λ, α are positive definite parameters signifying the strength of perturbation and type of singularity, has been of considerable interest for over three decades. This has a number of interesting physical as well as mathematical properties and found applications in chemical, nuclear and particle physics. For example, this represents a prototype of the so-called Klauder phenomenon [1, 2], i.e., even after the perturbation is completely turned off (λ → 0), some interaction still remains. The Rayleigh-Schrödinger perturbation series diverges following the relation n ≥ 1/(α − 2), with n denoting the order of perturbation. Another interesting feature is that in the region of α ≥ 5/2, this gives rise to supersingularity, i.e., the potential is so singular that any nontrivial correction to the energy (matrix element) diverges [1]. A fascinating property from a purely mathematical point of view is that there is no dominance of either of the two terms r 2 and λ/r α for the extreme values of λ. In other words, there is a tug-of-war between the two terms. Typical shape of the potential for a few parameter sets are depicted in Fig. 1 Like most of the practical physical systems, exact analytical solution of the corresponding Schrödinger equation (SE) for arbitrary parameters for a general state (both ground and excited) is yet to be known and consequently approximations such as variational or perturbation methods must be used. Thus a large number of attempts have been made towards the exact and approximate calculation of eigenvalues and eigenfunctions employing analytic, seminumerical or purely computational methods [3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18]. As the Rayleigh-Schrödinger perturbation series for the eigenvalues of the operator H diverges [3], a modified perturbation theory to finite order was developed which gave good result for the ground-state eigenvalues of small λ. For large positive values of λ, approximate estimates were made through the large coupling perturbative expansion [3]. A resummation technique for weak coupling of the nonsingular SHO (α < 5/2) [4] was developed. Exact and approximate solutions are also given for the lowest states of a given angular momentum [6]. A WKB treatment has be...