So, I wrote a generic nbody ODE solver, in pure python. It’s likely quite slow, but it works – seemed like too much work at the time to interface to scipy’s ODE solvers, which seemed designed for smaller, less coupled systems.
I may have misjudged, though which is why I’m noting this link http://www.scipy.org/Developer_Zone/Recipes/ODE_Integration_on_Banded_Systems ,which has an example of what seems to be a large dimensional ODE system. Thing is, we still need to do things like update adjacency lists and distances, and recompute summation densities between steps. Worth a look, though.