We present a method based on interpolating moving least-squares (IMLS) that is designed for efficient and accurate local fitting of discrete energy values to provide global representations of potential energy surfaces (PESs) for many-atom systems. We have demonstrated the method with one-, two-, and three-dimensional fits of the HN2 --> N-2 + H PES. To allow for extensive fitting and testing, the analytical PES developed by Koizumi et al. [Koizumi, H.; Schatz, G. C.; Walch, S. P. J. Chem. Phys. 1991, 95, 4130] was used to generate energy values. Unlike the modified Shepard method this approach does not require derivatives, thus it can be used to fit energies computed using highest-level quantum chemistry methods for which forces are not directly obtainable. This fitting scheme accurately describes the PES, is not computationally time-consuming, can be improved by using higher degrees and larger numbers of basis functions, and is straightforward to apply. Extension to many-dimensional PESs seems promising.