||A general method for iteratively fitting the coefficients of a Taylor-series expansion of the potential-energy surface for a polyatomic molecule to the observed transition frequencies and rotational constants is presented. This approach utilizes the efficiency of fourth-order Van Vleck perturbation theory for calculating these properties, scaled to the results of converged variational calculations. Three fits to the transition frequencies and rotational constants obtained from absorption and stimulated emission pumping spectra of HCN are presented. Comparison of these potentials sheds light on the uniqueness of fit potential surfaces, given this set of observables. We fit the vibrational energies for 71 states with a mean absolute deviation of 0.69 cm-1.