$\renewcommand{\Re}{\operatorname{Re}}$ $\renewcommand{\Im}{\operatorname{Im}}$ $\newcommand{\bR}{\mathbb{R}}$ $\newcommand{\bC}{\mathbb{C}}$ $\newcommand{\bZ}{\mathbb{Z}}$ $\newcommand{\const}{\operatorname{const}}$ $\newcommand{\sgn}{\operatorname{sgn}}$ $\newcommand{\rank}{\operatorname{rank}}$
It is known (see f.e. S. Chandrasekhar, Ellipsoidal Figures of Equilibrium that gravitational potential of ellipsoid $\mathcal{E}:=\bigl\{ (x,y,z)\colon \frac{x^2}{a^2}+\frac{y^2}{b^2}+\frac{z^2}{c^2}<1\bigr\}$, with the uniform density $\rho$ is \begin{gather} V(x,y,z)= -\pi abc\rho \int_\mu^\infty \left( 1-\frac{x^2}{a^2+\lambda} -\frac{y^2}{b^2+\lambda}-\frac{z^2}{c^2+\lambda}\right)\frac{d\lambda}{\Delta} \label{eq-L4.B-1} \end{gather} with $\mu=0$ as $(x,y,z)\in \mathcal{E}$ and with $\mu$ the positive root of equation \begin{gather} \frac{x^2}{a^2+\mu} +\frac{y^2}{b^2+\mu}+\frac{z^2}{c^2+\mu}=1. \label{eq-L4.B-2} \end{gather} as $(x,y,z)\notin \mathcal{E}$ and with \begin{gather} \Delta=\sqrt{(a^2+\lambda)(b^2+\lambda)(c^2+\lambda)} \label{eq-L4.B-3} \end{gather} Note that $\pi abc\rho =\frac{3}{4}M$ where $M$ is the mass of $\mathcal{E}$.
We are interested in the case of external point $(x,y,z)\notin \mathcal{E}$ when \begin{gather} a=b, \qquad c<a, \quad \varepsilon:=\frac{a^2}{c^2}-1 \ll 1. \label{eq-L4.B-4} \end{gather} Then \begin{gather} V(x,y,z)= \frac{3}{4}M \bigl(A + B\varepsilon + O(\varepsilon^2)\bigr) \label{eq-L4.B-5} \end{gather} where $A$ is the same integral with $a$ replaced by $c$ everywhere and $\mu$ replaced by $\mu^*=x^2+y^2+z^2-c^2$ (and one can prove easily that $A=\frac{4}{3} (x^2+y^2+z^2)^{-\frac{1}{2}}$) and
\begin{multline} B= -\int_{\mu^*}^\infty \frac{c^2(x^2+y^2)}{(c^2+\lambda)^{\frac{7}{2}}} d\lambda + c^2\int_{\mu^*}^\infty \left(1-\frac{x^2+y^2+z^2}{c^2+\lambda}\right)\frac{d\lambda}{(c^2+\lambda)^{\frac{5}{2}}}\\ +\left(1-\frac{x^2+y^2+z^2}{c^2+\lambda}\right) \Bigr|_{\lambda=\mu^*}\varepsilon^{-1}(\mu-\mu^*) \end{multline} with the last term vanishing.
Then $B=B_1+B_2$ with \begin{align*} B_1=-&\frac{2}{5}c^2(x^2+y^2) (c^2+\mu^*)^{-\frac{5}{2}}= -\frac{2}{5} c^2(x^2+y^2)(x^2+y^2+z^2)^{-\frac{5}{2}},\\ B_2=-&c^2\frac{2}{3}(c^2+\mu^*)^{-\frac{3}{2}}+\frac{2}{5}c^2(x^2+y^2+z^2)(c^2+\mu^*)^{-\frac{5}{2}}=-\frac{4}{15}c^2(x^2+y^2+z^2)^{-\frac{3}{2}} \end{align*}
Therefore \begin{multline} V(x,y,z)= - M (x^2+y^2+z^2)^{-\frac{1}{2}} - \\ M\Bigl(\frac{3}{10} c^2(x^2+y^2)(x^2+y^2+z^2)^{-\frac{5}{2}}+\frac{1}{5}c^2(x^2+y^2+z^2)^{-\frac{3}{2}}\Bigr)\varepsilon + O(M\varepsilon^2). \label{eq-L4.B-6} \end{multline}
$a=6,378.137\, km$, $c=6,355.752\, km$, $\varepsilon\approx 1/150$ rather small but not very-very small.
Therefore Low Eearth Orbits (few hundred km over Earth) of artificial satellites are indeed affected by the shape of Earth. Higher orbits, in particular, Geostationary Orbits are much less affected, and Moon is much–much less affected. Indeed, for them perturbation of $V=-MR^{-1}$ is $O(M\varepsilon R^{-3})$ where $R$ is the radius of orbit and therefore $\varepsilon$ is replaced by $\varepsilon c^2R^{-2}$.
Sure in reality $\rho$ is not uniform and Earth is not exactly ellipsoid of revolution but those are rather minor corrections.