  (Abridged:) Schwarzschild's method was used to construct equilibrium
solutions to the collisionless Boltzmann equation corresponding to a Plummer
sphere, which were compared with analytical results to test the robustness of
the numerical method and its efficiency in probing the degeneracy of the
solution space. The method was then used to construct triaxial equilibria with
no known analytical solutions and to study their nonuniqueness. The model that
was studied, the triaxial Dehnen potential, contains a central density cusp and
admits both regular and chaotic orbits. It was found that, for a weak density
cusp, self-consistent models do not exist if the chaotic orbits are assumed to
be completely mixed so as to yield a time-independent building-block: only the
innermost 65% of the mass can be mixed. In these inner regions, it is possible
to obtain alternative solutions that contain significantly different numbers of
chaotic orbits, yet yield (at least approximately) the same mass density
distribution. However, it is not clear whether these solutions are truly
time-independent, since the unmixed chaotic orbits in the outer regions, which
do not sample an invariant measure, can cause secular evolution. When using
Schwarzschild's method, one must be very careful to sample accessible phase
space as comprehensively and as densely as possible, while ensuring that each
orbit is a truly time-independent building block. Some of the numerical
equilibria were sampled to generate initial conditions for N-body simulations
to test the structural stability of the models. Preliminary work showed that no
catastrophic evolution takes place, but there is a weak tendency for the
configuration to become more nearly axisymmetric over several dynamical times.
It is not yet clear whether this tendency is real or a numerical artifact.
