Start by writing down the Schwarzschid metric for a spherically symmetric, nonrotating body:
<br />
ds^2 = -\left( 1 - \frac{2GM}{r} \right) \textrm{d} t^2 + \left( 1 - \frac{2GM}{r} \right)^{-1} \textrm{d} r^2 + r^2 ( \textrm{d} \theta^2 + \sin^2(\theta) \textrm{d} \varphi^2 ) \textrm{.}<br />
To get the geodesic equations from this, you can either calculate all the Christoffel symbols (long and tedious, but straightforward) or use the calculus of variations (my favorite method). In either case, I'd suggest setting e^{f(r)} = \left( 1 - \frac{2GM}{r} \right) for ease of notation. You'll end up with four very ugly-looking equations (one for each of t, r, \theta, and \varphi). By symmetry, you can set \theta = \frac{\pi}{2} (i.e., \theta' = 0). The t equation will then give you e^{f(r)} t' = E, where E is a constant of motion. The \varphi equation will give r^2 \varphi' = L, where L is another constant of motion ("angular momentum"). The r equation gives
<br />
e^{-f(r)} (r')^2 + \frac{L^2}{r^2} - e^{-f(r)} E^2 = \epsilon \textrm{,}<br />
where \epsilon is another consant of motion. For massive particles, you can set \epsilon = -1 without loss of generality (why?). Then you can rearrange the above equation to get
<br />
\frac{1}{2} (r')^2 + V(r) = \mathcal{E} \textrm{,}<br />
where V(r) = \frac{1}{2} \left(1 + \frac{L^2}{r^2} \right) e^{f(r)} is a "potential energy" and \mathcal{E} = \frac{1}{2} E^2 is a "total energy." This looks exactly like the form of a classical equation of motion, and you can analyze it using similar techniques. In particular, determining values for E, L, and r_{0} that lead to bounded or unbounded paths is straightforward.