I'm simulating the motion of a body in a Geocentric orbit by integrating the equation of motion.
I'm using the following equations for the acceleration of an object in a gravitational field defined by the central mass given by $\mu$ and an axisymmetric oblateness described by $J_2$:
Motion will be a Keplerian orbit, perturbed by the effects of $J_2$.
For increased accuracy, I should also consider other effects of Earth geometry, i.e J22, J3, J4. How could I modify the equation?
I would appreciate if you give links to papers.

+1While I've touched on the spherical harmonics of the Geopotential in this answer this new question needs a more thorough and mathematical answer than I can easily provide. Hopefully a "gravity person" will be able to address this better. – uhoh Feb 12 '18 at 08:23