In order to investigate the stability of stellar systems, we discuss a numerical method that uses an analytic distribution function to describe a stellar system in equilibrium and "perturbation particles" to represent departures from that equilibrium state. Thus, all the particles are used only to represent the perturbation, and statistical fluctuations due to the finite number of particles are much less severe than in full N-body codes. We provide a general description of the method, recipes for particular aspects of its implementation, and an example of its application to a simple model with known analytical solution.