A numerical method for studying collisionless stellar systems has been examined. The method exploits Liouville's theorem directly: the trajectories of phase particles are traced back to their initial positions to determine the evolution of the distribution function in phase space. Numerical results are presented, and the advantages and disadvantages of the method are discussed.