I want to produce some realistic figures showing the spatial trajectories of test particles in a Schwarzschild spacetime. For instance, I'd like to start a massive test particle at aponegricon (how often do you get to use that word!?) in an orbit that Kepler and Newton would have predicted to be elliptical, but that in fact crosses the event horizon. Basically if I can produce a list of (r,phi) pairs, I can produce such a plot. The strategy I tried initially was to use the equation of motion expressed in terms of r and phi, as described here: http://en.wikipedia.org/wiki/Schwarzschild_geodesics#Orbits_of_test_particles . We have an equation of the form (dr/dphi)^2=f(r), which I tried to integrate numerically. However, numerical integration of this equation is kind of a hassle. Since it gives you the square of dr/dphi, you have to decide on the sign, and the code has to figure out when you've reached a turning point and its time to flip the sign. In the case of a circular or near-circular orbit, I had an especially hard time making my code smart enough to get this right. Also, because of numerical errors you can have f(r)<0, which you then have to deal with, and you have to distinguish this from the case where the solution actually doesn't continue to larger phi, e.g., a hyperbolic orbit. The WP article also describes an exact solution using elliptic functions. This is a possible approach, although I would have to figure out if I can find an implementation of these functions in a programming language that's convenient. (I don't own a copy of Mathematica, but I do use Maxima sometimes.) A third approach would simply be to use (t,x,y) coordinates and apply the geodesic equation. I would think that this would be more numerically well-behaved than the (r,phi) version. Does anyone have any suggestions as to the easiest way to go about this?