The simplest and most efficient way of generating a set of N particles following some a priori density distribution (or corresponding gravitational potential) is to make use of cumulative functions.

Imagine you wish to generate positions for N particles following some arbitrary spherical density profile rho(r). The density does only depend on the radius, but not on the other 2 spherical coordinates Phi and Theta. So because of this simple geometry, you already know that each particle, for a fixed radius, has a uniform probability to be at an angle Phi, and Theta (in spherical coordinates). This means that if you already know the radius of each particle, you could then easily distribute them using spherical coordinates by randomly choosing a number between 0 and 2*Pi for Phi, and again one draw for Theta. You are then left with chosing the radius of each particle. This is where the cumulative function enters.

You can easily calculate the integrated mass M(r) for the model you wish to generate. M(r) is the mass within the sphere of radius r. So that if your model has a finite mass Mt (and only when this condition is satisfied), Mn(r) = M(r) / Mt is the normalized integrated mass, which is 0 for r=0, and 1 for r=infinity. Then one particle has a probability Mn(r) to be between a radius of 0 and r. You can then easily reverse the argument and generate N radius for N particles so that they follow the right density profile: you randomly choose a number X between 0 and 1, and set up the radius r such as Mn(r) = X. You then associate the corresponding radius r to that particle. This guarantees that after generating N particles, the integrated mass of your model will be consistent with M(r), so that the density profile of your model will naturally be consistent with rho(r).

The process to generate the positions for N particles is then easily described:
 you first derive Mn(r) for your favorite rho(r)
 you then randomly draw numbers X_i between 0 and 1 for your N particles (using a uniform distribution), where i goes from 1 to N.
 you associate the corresponding radii to each of them (Mn(r_i) = X_i).
 you then draw N times a random number between 0 and 2*Pi for Phi, and the same for Theta.
 you end up with r, Theta, Phi for your N particles, consistently with the initially given density distribution rho(r).

This illustrates how to generate a set of N bodies for a spherical distribution. This can obviously be translated for any other distribution where the density distribution has a form which allow coordinates to be independent from each others. For example: you can easily generate N particles in a Toomre disk by first assuming a two-dimensional distribution (zero thickness), and using the cumulative function for such a radial profile, and then thicken the disk using some random distribution in along the vertical axis (uniform, gaussian, etc...). You should also have a look at the detailed explanation provided in the Appendix of the original paper by Aarseth, Henon, Wielen 1974, A&A 37, 183, or in the dialog set up by Piet Hut and Jun Makino (Kali).

Many examples exist to illustrate this, but the best is to have a look at codes available on the web. You can for instance have a look at the software distributed by Joshua Barnes (e.g. The Zeno programs at Zeno, or his latest tree code at Treecode).