The contour integral method of Hunter & Qian is applied to axisymmetric galaxy models in which the distribution function (DF) is of the form f=f(E,L_z), where E and L_z are the classical integrals of motion in an axisymmetric potential. A practical way to construct the unique even part of the DF for such systems is presented. It is applied to models, both oblate and prolate, in which the mass density is stratified on similar concentric spheroids. The spheroids with scale-free densities are discussed in detail. These provide useful approximations to the behaviour of more realistic models in the limit of small and large radii. The self-consistent case is treated, as well as the case in which there are additional contributions to the potential from a central black hole or dark halo. The range of axis ratios and density profile slopes is determined for which spheroidal power--law cusps with a central black hole have a physical two--integral DF. More generally, the f(E,L_z) DFs are discussed for a set of spheroidal `(alpha,beta)-models', characterized by a power-law density cusp with slope alpha at small radii, and a power-law density fall-off with slope [alpha + 2 beta] at large radii. As an application, the DF is constructed for the (alpha,beta) model with a 1.8 x 10^6 solar mass black hole used by van der Marel et al. to interpret their high spatial resolution spectroscopic data for M32. The results confirm that the model fits the observed line-of-sight velocity profiles remarkably well. The model is used to calculate predictions for future spectroscopic observations with the HST.