I am trying to derive the rectangular components of acceleration for a satellite in orbit, with Earth oblateness in consideration, in order to use the RK4 method to find the updated position and velocity of the satellite. I am aware that there are equations already out there, but I want to know how they were derived. However, my solutions don't exactly look like those equations - more importantly, the power of r in those sets of equations and in the wikipedia force equations are not the same. These equations have $\frac{1}{r^7}$, while mine have $\frac{1}{r^6}$
I will show how I derived the equations, please let me know if I've made an error!
I am using the following equation that I grabbed from the NASA Flight Test Cases document and Fonte 1993 (Implementation of 50x50 Gravity Field Model):
$$U = -\frac{\mu }{r}\left [ 1+\sum_{n=2}^{\infty}\sum_{m=0}^{n}\left(\frac{ R_{e}}{r}\right)^{n}P_{n,m}(sin\phi )\left(C_{n,m}cos(m\lambda) +S_{n,m}sin(m\lambda ) \right)\right]$$
I knew that for zonal harmonics $m=0$, so the equation looks like $$U = -\frac{\mu }{r}\left[1+\sum_{n=2}^{\infty}\left(\frac{R_{e}}{r}\right)^{n}P_{n,0}(sin\phi )(C_{n,0}cos(0) )\right]$$
I will then solve for $n=2$ knowing that $C_{n,0}=J_{2}$
There is a section below equation 10 under the deviations of Earth's gravitational field from that of a homogenous sphere section of the Geopotential Wiki page that says
$$u=\frac{J_{2}P_{2}^{0}(sin\theta )}{r^3}=J_{2}\frac{1}{r^3}\frac{1}{2}(3sin^{2}\theta-1)=J_{2}\frac{1}{r^5}\frac{1}{2}(3z^{2}-r^{2})$$
So I deduced that
$$P_{2,0}=\frac{1}{r^2}\frac{1}{2}(3z^{2}-r^{2})$$
So the equation looks like
$$U = -\frac{\mu }{r}\left[1+\left(\frac{R_{e}}{r}\right)^{2}(J_{2} )\left(\frac{1}{r^2}\frac{1}{2}\left(3z^{2}-r^{2}\right)\right)\right]$$
I then multiplied everything out and put the equation in this form
$$U = -\mu\left[\frac{1 }{r}+\frac{3}{2}J_{2}R_{E}^{2}z^{2}\frac{1 }{r^5}-\frac{1}{2}J_{2}R_{E}^{2}\frac{1 }{r^3}\right]$$
Knowing that $F_{x} = -\frac{\partial U}{\partial x}$, and using $\frac{\partial }{\partial x}r^{-5}=\frac{\partial }{\partial x}(x^{2}+y^{2}+z^{2})^{-5}=-\frac{10x}{r^6}$ etc., and then rearranging, I got the following for the x (and y) component:
$$F_{x} = -\frac{\partial U}{\partial x} = -\frac{\mu x}{r^2}\left[2+\frac{3J_{2}R_{E}^{2}}{r^2}\left(\frac{5z^{2}}{r^2}-1\right)\right]$$
Which looks somewhat familiar to this, but is not it exactly. I've tried doing that math three times now and get the same answer.
Sidenote: I asked OP of that thread, and he said that he got those equations from Section 7A of this competition
Can anyone tell me what I'm doing wrong?