We discuss a numerical method for investigating the stability of stellar systems that uses an analytic distribution function to describe a stellar system in equilibrium and perturbation particles to represent departures from that equilibrium state. The particles are used only to represent the perturbation, and statistical fluctuations due to the finite number of particles are therefore 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.