I'm trying to model numerically the effects of the Earth's oblateness on satellite orbits. I'm seeing precession of the orbit, but not in the plane that I'm expecting. Here's my corrected grav. potential energy: $$U=-\frac{GMm}{r}+\frac{GMma^2}{2r^3}(2cos^2\theta -1)$$
where $\theta$ is the colatitude (measured from the pole, not the equator).
I think that the problem may be in getting from $U$ to $F$ (which is what I really need for the simulation). I'm differentiating w.r.t. $r$ - do I also need to do something about the $\theta$ variation when I go from $U$ to $F$ here?