In this work we develop and test a methodology for the generation of Born-Oppenheimer potential energy surfaces (PES) for use in vibrational structure calculations. The method relies on the widely used restricted-mode-coupling expansion of the fully coupled potential surface where only up to n or less vibrational coordinates are coupled in the potential. Low-order derivatives of the energy are then used to extrapolate the higher mode-coupling potential terms; derivative information is thus used in a convenient way for the evaluation of higher mode couplings avoiding their explicit calculation on multidimensional grids. The formulation, which is a variant of the popular modified Shepard interpolation, is general for any extrapolation of (n+p)-mode-coupling terms from n-mode couplings and can be applied to the energy or any other molecular property surface for which derivative information is available. The method depends only on analytical parameter-free weight functions that satisfy important limiting conditions and control the contribution from each direction of extrapolation. The procedure has been applied on a representative set of 13 molecules, and its accuracy has been tested using only gradients and using both gradients and Hessians. The results provide evidence for the importance of higher mode couplings and illustrate the cost efficiency of the proposed approach.