We proposed a global spin-dependent potential energy function which fits the description of the non-relativistic quark model. The spin dependencies of our model result from the harmonic approximation of the Gaussian function in the spin-spin component. The approximation enables us to obtain analytical solutions of the bound states for the heavy mesons singlet and triplet states using the conventional Nikiforov-Uvarov method. The obtained mass spectra and expectation values for the low through high-lying quantum states are in agreement with the results obtained for other potential functions in the existing literature. The splitting of the masses occurs at a distance greater than the obtained characteristic radius of the heavy mesons which were found to be within the experimentally observed 1fm. The proposed potential model may be extended to obtain other meson spectroscopic parameters such as the decay constants, wave function at the origin and electromagnetic radiative transitions.