The nature of the integrals of motion in Staeckel potentials, for which the equations of motion separate in ellipsoidal coordinates, is elucidated. The problem of fitting a general potential with one of Staeckel form is considered, both locally and globally. In local fitting the potential is expanded around an equilibrium point and this expansion is compared term by term with the similar expansion of a Staeckel potential. This is done explicitly for potentials with three reflection symmetries. Expansions are given for the integrals of motion in the best-fitting Staeckel potential. They may be expressed directly in terms of the expansion coefficients of the potential that is fitted. The results of local fitting are applied to the gravitational potentials of ellipsoidal density distributions, and also to ellipsoidal potentials. It is shown that there is a unique inhomogeneous density distribution stratified on similar ellipsoids of arbitrary axis ratios with a potential that is exactly of Staeckel form. A method is presented for the global fitting of a triaxial potential with one of Staeckel form. The foci of the ellipsoidal coordinates play an important role in this procedure.