Finding the Right Initial Conditions for a Three-Body Problem Periodic Solution The planar restricted circular three-body problem is about computing the trajectory of a small mass, usually a speccraft, that is affected by two larger masses, Earth and the Moon. The differential equation for this problem is given as 
$$ 
y_1'' = y_1 + 2y_2' - \mu_2 \frac{y_1+\mu_1}{D_1} - \mu_1 \frac{y_1-\mu_2}{D_2}
$$
$$ 
y_2'' = y_2 - 2y_1' - \mu_2 \frac{y_2}{D_1} -\mu_1 \frac{y_2}{D_2}
$$
where
$$\mu_1 = \frac{m_1}{m_1+m_2}$$
$$\mu_2 = 1-\mu_1$$
$$ D_1 = ((y_1+\mu_1)^2 + y_2^2 )^{3/2}$$
$$ D_1 = ((y_1-\mu_2)^2 + y_2^2 )^{3/2}$$
For these equations for a certain initial conditions a so-called Arenstorf orbit can be computed which is periodic. The initial conditions are given as (for example in Solving Ordinary Differential Equations I, Nonstiff Problems, pg. 130),
$y_1(0) = 0.994$, $y_1'(0)=0$, $y_2'(0) = -2.0015851063790825224053786222$
Every textbook, lecture note repeats these values, which produces, after integration, a nice orbit around earth and moon. But noone seems to know where they come from. Were they computed analytically? Or a trial and error approach was used to find them, or a combination of both? 
If someone has code that can compute these values it would be much appreciated. 
Thanks,
 A: One thing to note is that, while you have a system of two equations, since each one is a second order equation, you have a fourth order system. That means that the phase space is 4D, so the separatrices are 3D manifolds in 4D space, and orbits are 1D manifolds (curves) in 4D space. You can project the orbits onto the plane, but they lose their property of being non-self-intersecting. Thus, most of the math that goes into this cannot be visualized, and so even the simplest coupled oscillators can provide a real challenge from a mathematical standpoint. Most of the theory on these things has been developed in the last 50 or so years, so keep that in mind when reading this. Also keep in mind that I'm a mathematician who studies dynamics, not a physicist or engineer, so while I have experience with numerical methods, my background is in analytical and theoretical methods.

For most strongly nonlinear systems, numerics are the only way (currently) to determine precisely the behavior, especially with respect to various parameters. For varying parameters, we usually refer to bifurcation theory. For a given set of parameters, we can identify the basins of attraction, which are the IC regimes where the system tends to each behavior. The boundaries between these regimes are called separatrices (singular: separatrix), or invariant manifolds. In autonomous systems, these manifolds are nice and simple (meaning not fractals) and each basin is a connected domain (two IC's in the same regime always have a curve connecting them entirely in the regime).
So, what we can do to reduce the trial and error to reasonable amounts is to try to find the fixed points and their linearized behavior, approximate the separatrices to get a partitioned space, then choose some sample IC's within each partition to try and determine the behavior of the system in that regime. You may stumble upon some points which generally move to fixed points, some which move to periodic orbits, some which move to quasiperiodic orbits, and some which move to chaos. Within the regime where a cycle happens, there must be some fixed point. You can transform that fixed point to the origin, then get the amplitude/phase system by assuming $y_i = r_i \cos(t+\psi_i)$ for $i=1,2$, then use variation of parameters and averaging to get an approximation of the cycle. 
For the linear stability analysis, see any undergraduate book on Linear and Nonlinear ODE's. Approximation methods for invariant manifolds usually come in the form of some power series methods, but you can search for those. Once you've determined some approximate partitioning, just pick a bunch of initial conditions, then numerically solve the equations (using, say, ode45 with MATLAB) for each of the IC's and plot the results on the $(y_1,y_2)$-plane. For the variation of parameters and averaging, you can look at some books on Nonlinear Oscillations or Perturbation Methods, for example by Ali Nayfeh.
A: A professor answered in private communication that "the solutions are found numerically by solving a two-point boundary value problem with periodic boundary conditions. The two-point boundary problem is solved using Newton's method. There is some trial and error involved in finding a suitable starting point for the Newton iterations". 
Based on this, I did a search, and found this paper. How to solve the boundary problem between point A and B? The author says guess an initial velocity $v_i$, and since starting position known, all initial value are known. Integrate, if the result is far from point B, not within alloted time, update $v_i$ according to a gradient formule, integrate again.
