Besides a variety of programs of the kind of mkplummer, mkommod, mkexpdisk etc., models can also be generated by calculating appropriate tables (containing a run of density, potential and radius) and feeding them into programs which translate such tables into a snapshot. An example of a program in the making is anisot.
An alternative way is to use a package such as Mathematica to integrate the differential equations. With the following example we leave this as an exercise to the reader:
<< NumericalMath/RungeKutta.m rpsipsiprime = RungeKutta[{psiprime, -(2/r)psiprime + Exp[-psi]}, {r, psi, psiprime}, {10^-5, 10^-10/6, 10^-5/3}, 50, 10^-8, MaximumStepSize->0.5]; rpsi = rpsipsiprime[[Table[n,{n,Length[rpsipsiprime]}],{1,2}]]; r = rpsi[[Table[n,{n,Length[rpsi]}],1]]; psi = rpsi[[Table[n,{n,Length[rpsi]}],2]]; rho = Exp[-psi]; rrhotranspose={r,rho}; rrho = Transpose[rrhotranspose]; logrlogrho = Map[Log, rrho, {2}]; c = Log[10.]; log10rlog10rho = logrlogrho / c; shortrrho = Take[rrho,30]; PlotODESolution[shortrrho, 1, 2, AxesLabel-> {"r", "rho"}, PlotLabel -> "Isothermal Sphere"] PlotODESolution[log10rlog10rho, 1, 2, AxesLabel-> {"log r", "log rho"}]