Consider the following setup:

Two massive bodies of mass $M$ are **fixed** at the positions $\vec{A}=(d, 0)$ and $\vec{B}=(-d, 0)$.
![setup][1]

Now imagine a test particle $p$ with initial position $\vec{r}_0=(x_0,y_0)$. It's coordinates will change under the influence of the gravitational potential. Its coordinates are $x(t)$ and $y(t)$. If the test particles trajectory goes through $A$ (or $B$) the initial position lies in the basin of attraction of the mass at point $A$ (or $B$).

**Question**: Does $x_0 \gt 0$ imply that $p$ will eventually come arbitrarily close to $A$ (**and not** $B$)?

**Or simpler yet**: Does $x_0 \gt 0$ imply that $x(t)\gt0$ for all $t$?


---
Background:

I'm currently working on a (2 dimensional) simulation of gravity which visualizes gravitational basins of attraction.

The simulation works as follows:

$n$ massive bodies are given an intial position. They are treated as **static** and will never change position during the simulation. The bodies are given the **same mass**.

Next, the simulation space is filled with particles. Their initial velocity is set to zero. Their trajectories are calculated via numerical integration. 

When a collision between a particle and one of the massive bodies is detected, the initial position of the particle is assumed to lie in the basin of attraction of that body (with which it collided). A collision occurs when the distance to a body is below a certain threshold.

---

Simulating the above described setup gave the following picture:

[![basins for two bodies][2]][2]

Where blue pixels represent initial positions in the basin of attraction of the mass at point $A$, red pixels those in the basin of the mass at point $B$ and black pixels in neither. (A thin strip of black pixels lie on the line $x=0$, where they oscillate around the origin)

---

EDIT:

I tried solving the problem analytically using Lagrange's equations of the second kind. I set:
$$T=\frac{m}{2}(\dot{x}^2+\dot{y}^2)$$
$$V=-GM\left[\frac{1}{\sqrt{(x-d)^2+y^2}}+\frac{1}{\sqrt{(x+d)^2+y^2}}\right]$$
And obtained the following equations of motion:
$$\ddot{x} = -\frac{GM}{m}\left[\frac{x-d}{\left((x-d)^2+y^2\right)^{\frac{3}{2}}} + \frac{x+d}{\left((x+d)^2+y^2\right)^{\frac{3}{2}}}\right]$$
$$\ddot{y} = -\frac{GM}{m}\left[\frac{y}{\left((x-d)^2+y^2\right)^{\frac{3}{2}}} + \frac{y}{\left((x+d)^2+y^2\right)^{\frac{3}{2}}}\right]$$

Where $x$ and $y$ are the coordinates of the particle, $m$ the mass of the particle, $M$ the mass of the massive bodies and $d$ the distance (in the $x$-direction) from the origin of the massive bodies.

I converted this into a system of four first-order differential equations and tried solving it using Maple's *dsolve*, but it's been stuck *Evaluating* for two hours now, so I don't think it's going to finish...


  [1]: https://i.sstatic.net/Jx9kjl.png
  [2]: https://i.sstatic.net/IXXoJ.png