In predictive ab initio spectroscopy, the use of high-level theory and large basis sets makes each energy calculation computationally expensive. As a result, the energy points are not computed on a fine mesh. Continuous Born-Oppenheimer potential-energy functions are generated by a smoothing procedure. For example, cubic splines are used in the VIBROT/MOLCAS program and high-order polynomials are used in the MOLPRO program. Several smoothing options are available in the popular LEVEL program, by the late R. J. LeRoy. Different smoothing procedures produce different spectroscopic results. Even within a single model, such as splines, interpolation error may be significant for vibrational eigenvalues. The choice of grid points may also affect the results. For sparse data, which is a common situation, the differences between smoothing methods can be surprisingly large, as shown by the examples below. Here we propose a strategy to reduce this variability.