  We describe a practical method for constructing axisymmetric two-integral
galaxy models (with distribution functions of the form f(E,L_z), in which E is
the orbital energy, and L_z is the vertical component of the angular momentum),
based on Schwarzschild's orbit superposition method. Other f(E,L_z)-methods are
mostly based on solving the Jeans equations or on finding the distribution
function directly from the density, which often places restrictions on the
shape of the galaxy. Here, no assumptions are made and any axisymmetric density
distribution is possible. The observables are calculated (semi-)analytically,
so that our method is faster than most previous, fully numerical
implementations. Various aspects are tested extensively, the results of which
apply directly to three-integral Schwarzschild methods. We show that a given
distribution function can be reproduced with high accuracy and investigate the
behaviour of the parameter that is used to measure the goodness-of-fit.
Furthermore, we show that the method correctly identifies the range of cusp
clopes for which axisymmetric two-integral models with a central black hole do
not exist.
