Correcting for oblateness of the Earth [duplicate]

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?

Yes, the gradient in spherical coordinates contains angular terms. It has to, because the force need not point to the center. The relevant equation is $$\Delta U=\frac {\partial U}{\partial r}e_r+\frac {\partial U}{r\partial \theta}e_\theta+\frac {\partial U}{ r \sin \theta \partial\phi}e_\phi$$ You have no $\phi$ variation, so can ignore the last term.