We present a novel, least-squares fitting approach to obtain a representation of a potential energy surface using potential, gradient, and Hessian data. The method is described in detail and then tested for H2O in two ways. In the first test a global, analytical potential is used to generate the data at 7 and 11 configurations. A comparison of the accuracy of the fit against the exact surface is made, as is a comparison of low-lying vibrational states. In the second test, Density Functional theory (DFT) calculations of the potential, gradient, and Hessian are performed at 7 and 11 configurations to obtain fits. The predictions of the fits are compared to 125 new DFT calculations of the energies and a conventional fit to them, both directly and in vibrational calculations.