  We used the N-body code of Hernquist and Ostriker (1992) to build a dozen
cuspy ({\gamma}\approx 1) triaxial models of stellar systems through
dissipationless collapses of initially spherical distributions of 10^6
particles. We chose four sets of initial conditions that resulted in models
morphologically resembling E2, E3, E4 and E5 galaxies, respectively. Within
each set, three different seed numbers were selected for the random number
generator used to create the initial conditions, so that the three models of
each set are statistically equivalent. We checked the stability of our models
using the values of their central densities and of their moments of inertia,
which turned out to be very constant indeed. The changes of those values were
all less than 3 per cent over one Hubble time and, moreover, we show that the
most likely cause of those changes are relaxation effects in the numerical
code. We computed the six Lyapunov exponents of nearly 5,000 orbits in each
model in order to recognize regular, partially and fully chaotic orbits. All
the models turned out to be highly chaotic, with less than 25 per cent of their
orbits being regular. We conclude that it is quite possible to obtain cuspy
triaxial stellar models that contain large fractions of chaotic orbits and are
highly stable. The difficulty to build such models with the method of
Schwarzschild (1979) should be attributed to the method itself and not to
physical causes.
