An adapted formulation of the grid-based Numerov approach to solve Schrödinger’s equation has been extended to four-dimensional systems and applied to problems in vibrational spectroscopy. The advantage of this novel implementation is that no approximations or assumptions with respect to the wavefunctions are required and the achieved accuracy at a given stencil size is only limited by the selected grid-spacing. To validate this approach the vibrational eigenstates and the corresponding wavenumbers of the linear molecules CO2, BeH2 and HCN were investigated. Based on a potential energy grid at CCSD(T)/cc-pVnZ (n = T, Q) level all fundamental wavenumbers could be reproduced with deviations of less than 1% of the experimental values. Furthermore, the hierarchical application of a one- to four-dimensional Numerov treatment provides detailed information about the increasing influence of inter-mode coupling in addition to the inherent account for anharmonic contributions. The associated vibrational states are eigenfunctions of the Numerov matrix eigenvalue problem. Since no vibrational basis functions have to be applied in this approach, the respective eigenstates are characterized using isosurfaces of the wavefunctions in the hyperplane x4=0. This enables the direct examination of the individual wavefunctions via suitable visualization tools and the assignment of the respective vibrational quantum numbers.