You have four big problems and two small problems. The big problems are that you are initializing the Earth's and Moon's initial position and velocity incorrectly. The initial distance between the Earth and Moon is off by a factor of 1.0123, as is the initial relative velocity. The small problems are (1) an incorrect value for the Earth-Moon semi-major axis and (2) your use of $G$, $m_{moon}$, and $m_{earth}$.
The big problems were the key cause of your larger than expected apogee distance.
The small problems: You should want to fix those as well.
Issue #1: You are using 384,399 km as the length of the semi-major axis of the the orbit of the Moon. That is incorrect. That value is the inverse sine parallax of the Moon, the inverse of the mean value of the inverse of the distance. A better value for the semi-major axis length is 385,000 km, which is the mean distance between the Earth and Moon. An even better value is 384,748 km, from Chapront-Touzé, M., & Chapront, J. (1983). The lunar ephemeris ELP 2000. Astronomy and Astrophysics, 124, 50-62.
Issue #2: You are using the product $G (m_{earth}+m_{moon})$. You have obliterated your accuracy when you do that. Solar system astronomers instead use what are called the "standard gravitational parameters" to describe the masses of the Sun, the planets, and our Moon. Conceptually, this is just the product $\mu_{body} = G m_{body}$. There's a big difference, though, between using $\mu_{body}$ and $G m_{body}$. Scientists know many of those gravitational parameters to six places of accuracy or more. On the other hand, the gravitational constant G is known a paltry four places. Even more importantly, when you use G and mass you are almost certainly using values that are inconsistent with one another. Use the standard gravitational parameters. You can find a list of them at this wikipedia article.
So what should you do?
Denoting
- $r_p=362600\,\text{km}$ as the distance between the Earth and Moon at perigee,
- $a= 384748\,\text{km}$ as the semi-major axis of the Earth and Moon about one another,
- $u_e=398600.4418\,\text{km}^3/\text{s}^2$ as the Earth's standard gravitational parameter,
- $u_m=4902.8000\,\text{km}^3/\text{s}^2$ as the Moon's standard gravitational parameter, and
- $v_p = \sqrt{(\mu_e+\mu_m)\left(\frac 2 {r_p} - \frac 1 a\right)}\,$ as the relative velocity between the Earth and Moon at perigee, per the vis-viva equation.
You need to emplace the Earth and Moon such that the distance and velocity between them are $r_p$ and $v_p$. Since you want the barycenter to be at the origin, one way to do this is
$$\begin{aligned}
r_{moon} &= \phantom{-} \frac {r_p} {1+\mu_m/\mu_e} \hat x \\
v_{moon} &= \phantom{-} \frac {v_p} {1+\mu_m/\mu_e} \hat y \\
r_{earth} &= -\frac {r_p} {1+\mu_e/\mu_m} \hat x \\
v_{earth} &= -\frac {v_p} {1+\mu_e/\mu_m} \hat y
\end{aligned}$$
Finally, you should compute your accelerations using the standard gravitational parameters rather than using G*M.
We can the path of any moon as an epicycloid: we treat the planet as though it moves in uniform circular motion with radius $r_p$, and the moon as though it moves in uniform circular motion about the planet with radius $r_m$. The net acceleration of the moon will be the vector sum of the planet's acceleration towards the sun (with magnitude $r_p \omega_p^2$) and the moon's acceleration towards the planet (with magnitude $r_m \omega_m^2$). Here, $\omega_p$ is the angular velocity of the planet about the sun, and $\omega_m$ is the angular velocity of the moon about the planet.
Since the planet's acceleration always points towards the sun, the only way for the moon to have a net acceleration vector pointing away from the sun is for its acceleration vector towards the planet to be of greater magnitude:
$$
r_m \omega_m^2 > r_p \omega_p^2,
$$
or in terms of the respective periods $T_p$ and $T_m$,
$$
\frac{r_m}{T_m^2} > \frac{r_p}{T_p^2}.
$$
Under these approximations, any moon for which this is satisfied will have some point of "concavity" in its orbit when it is between the Sun and its parent planet.
I whipped up some code using Mathematica and Wolfram's curated data to see which moons satisfy this condition. My initial findings are that only a few of the outermost moons of Jupiter, Saturn, and Neptune do not have points of concavity in their orbits (along with Earth's moon, of course.) This makes a certain amount of sense, since orbital periods of moons aren't significantly longer in the outer solar system than in the inner solar system, while the orbital periods of the planets definitely are.
However, I'm not sure that I've done the coding correctly, since I expected my code to return all moons of a given planet outside a certain radius and this isn't what happened. I will try to fix this and update once I have a better grasp on the real answer.
Best Answer
By the sounds of it you have made a mistake with the units. In fact, you should not be using SI units at all in your simulation; astronomical values in SI units vary by such huge orders of magnitude that they are often a source of floating point errors that can destroy trajectories.
You should instead use the astronomical system of units. Specifically, express your masses in solar masses, your lengths in the astronomical unit, and time in the mean solar day.
Your value of $G$ will then be the square of the Gaussian gravitational constant, i.e.
$$ G=k^2=0.0002959122083\,\mathrm{AU}^3\mathrm{D}^{-2}\mathrm{M}_\odot^{-1} $$
You can get position and velocity data for the planets, their moons, comets, and hundreds of thousands of asteroids, etc, from JPL HORIZONS. You need to connect to their servers via telnet and request the data.
Alternatively, as a starting point, if you know the distance $r$ and mass $m$ of a planet, then its orbital speed should be
$$ v = \sqrt{\frac{GM}{r}} $$
where $M=1\,\mathrm{M}_\odot$. If we ignore eccentricity then the direction should be tangential to the orbit (e.g. position it at $(r,0)$ and give it a velocity $(0,v)$).