There are lots of things you can do. It is a simple ODE, and can be solved straightforwardly using numerical methods in specific spacetimes. If you write the equation in Fermi-like coordinates based around the worldline of interest, it turns into something like a harmonic oscillator equation with variable (and possibly negative) spring constant.
More interestingly, all solutions are known analytically (at least in a finite region) in terms of the spacetime's world function (also known as Synge's function). This is a two-point scalar [itex]\sigma(x,x') = \sigma(x',x)[/itex] that returns one-half of the geodesic distance between its arguments. See Dixon, Proc. R. Soc. A 314, 499 (1970) or http://arxiv.org/abs/0807.1150" . Standard results in bitensor perturbation theory can be used to turn this solution into a power series near the initial point, if desired. Part of the exact solution also looks very nice in Riemann normal coordinates (without approximation).
Alternatively, one can write down an equation that leads to an iterative approximation scheme using integrals of the Riemann tensor. This is useful for cases where the Riemann tensor changes considerably but never gets particularly large. See the GR textbooks by Synge or de Felice and Clark.