We present a scheme for numerical simulations of collisionless
self-gravitating systems which directly integrates the Vlasov--Poisson
equations in six-dimensional phase space. By the results from a suite of
large-scale numerical simulations, we demonstrate that the present scheme can
simulate collisionless self-gravitating systems properly. The integration
scheme is based on the positive flux conservation method recently developed in
plasma physics. We test the accuracy of our code by performing several test
calculations including the stability of King spheres, the gravitational
instability and the Landau damping. We show that the mass and the energy are
accurately conserved for all the test cases we study. The results are in good
agreement with linear theory predictions and/or analytic solutions. The
distribution function keeps the property of positivity and remains
non-oscillatory. The largest simulations are run on 64^6 grids. The computation
speed scales well with the number of processors, and thus our code performs
efficiently on massively parallel supercomputers.Comment: 35 pages, 19 figures. Submitted to the Astrophysical Journa