The classification of short hydrogen bonds depends on several factors including the shape and energy spacing between the nuclear eigenstates of the hydrogen.Here, we describe the NuSol program in which three classes of algorithms were implemented to solve the 1D, 2D and 3D time independent nuclear Schrödinger equation. The Schrödinger equation was solved using the finite differences based Numerov’s method which was extended to higher dimensions, the more accurate pseudo-spectral Chebyshev collocation method and the sinc discrete variable representation by Colbert and Miller. NuSol can be applied to solve the Schrödinger equation for arbitrary analytical or numerical potentials with focus on nuclei bound by the potential of their molecular environment. We validated the methods against literature values for the 2D Henon–Heiles potential, the 3D linearly coupled sextic oscillators and applied them to study hydrogen bonding in the malonaldehyde derivate 4-cyano-2,2,6,6-tetramethyl-3,5-heptanedione.With NuSol, the extent of nuclear delocalization in a given molecular potential can directly be calculated without relying on linear reaction coordinates in 3D molecular space. Program summaryProgram title: NuSolCatalogue identifier: AEXO_v1_0Program summary URL:http://cpc.cs.qub.ac.uk/summaries/AEXO_v1_0.htmlProgram obtainable from: CPC Program Library, Queen’s University, Belfast, N. IrelandLicensing provisions: GPL 3.0No. of lines in distributed program, including test data, etc.: 195332No. of bytes in distributed program, including test data, etc.: 4489808Distribution format: tar.gzProgramming language: Python 2.7, C.Computer: PC.Operating system: Linux.RAM: 0.1–40 GBClassification: 16.1.External routines: FEAST v2.1 [1][2] (included in the distribution file).Nature of problem: Solving the 3D nuclear Schrödinger EquationSolution method: Grid solver based on Numerov’s method, Chebyshev collocation and sinc() DVR for 1D/2D/3D potential gridsRunning time: System dependent