Over the past decade, about 50 potential energy surfaces (PESs) for polyatomics with 4-11 atoms and for clusters have been calculated using the permutationally invariant polynomial method. This is a general, mainly linear least-squares method for precise mathematical fitting of tens of thousands of electronic energies for reactive and nonreactive systems. A brief tutorial of the methodology is given, including several recent improvements. Recent applications to the formic acid dimer (the current record holder in size for a reactive system), the H2-H2O complex, and four protonated water clusters [H+(H2O)n=2,3,4,6] are given. The last application also illustrates extension to large clusters using the many-body representation.