Solving Keplerian two-body problem numerically? I'm trying to solve the two-body problem numerically, setting up G, m1 and m2 =1.0.
 The masses are placed at the positions -10 and 10 respectively along the x-axis and gave them both 0 on the y-axis. I am having some real issue with the initial conditions fitting for a keplerian orbit with a given eccentricity (e=0.9)!
Also how would the initial condition change for different mass ratios?
 A: If we give the two masses initial velocities in the $\hat{y}$ and $-\hat{y}$ directions, perpendicular to their separation, we'll be starting at either periapsis or apoapsis.
The vis-viva equation for a binary system says that
$$v^2=G(m_1+m_2)\left(\frac{2}{r}-\frac{1}{a}\right)\tag1.$$
Here $r$ is the distance between the two masses, which is the magnitude of the separation vector
$$\mathbf{r}=\mathbf{r}_1-\mathbf{r}_2\tag2;$$
note that this is not the distance of either from their barycenter.
The speed $v$ is the magnitude of the relative velocity vector
$$\mathbf{v}=\frac{d\mathbf{r}}{dt}=\mathbf{v}_1-\mathbf{v}_2\tag3.$$
The semimajor axis $a$ is for the ellipse formed by $\mathbf{r}$, not the smaller ellipses formed by $\mathbf{r}_1$ or $\mathbf{r}_2$.
The ellipse formed by $\mathbf{r}$ has the form
$$r=\frac{a(1-e^2)}{1+e\cos\theta}\tag4$$
where $\theta=0$ at periapsis and $e$ is the eccentricity.
The value of $r$ at apoapsis is
$$r_a=a(1+e)\tag5.$$
So the relative speed at apoapsis is
$$\begin{align}
v_a&=\sqrt{G(m_1+m_2)\left(\frac{2}{r_a}-\frac{1}{a}\right)}\\
&=\sqrt{G(m_1+m_2)\left(\frac{2}{r_a}-\frac{1+e}{r_a}\right)}\\
&=\sqrt{\frac{G(m_1+m_2)}{r_a}(1-e)}\\
\end{align}\tag6.$$
Thus, given your initial data, we can find the initial relative speed.
To determine the initial speed of each mass, we use the fact that
$$\mathbf{r_1}=\frac{m_2}{m_1+m_2}\mathbf{r}\tag7$$
and
$$\mathbf{r_2}=-\frac{m_1}{m_1+m_2}\mathbf{r}\tag8$$
which follow from (2) plus the condition that the center of mass is at the origin:
$$\frac{m_1\mathbf{r}_1+m_2\mathbf{r}_2}{m_1+m_2}=0\tag9.$$
We thus have
$$\mathbf{v_1}=\frac{m_2}{m_1+m_2}\mathbf{v}\tag{10}$$
and
$$\mathbf{v_2}=-\frac{m_1}{m_1+m_2}\mathbf{v}\tag{11}$$
Putting in your numbers $G=m_1=m_2=1$, $r_a=(10)-(-10)=20$, and $e=9/10$, we find the initial speeds should be
$$v_1=v_2=\frac{1}{1+1}\sqrt{\frac{(1)(1+1)}{20}\left(1-\frac{9}{10}\right)}=\frac{1}{20}\tag{12}.$$
Trying this in Mathematica using

data = NBodySimulation[
"InverseSquare", {<|"Mass" -> 1, "Position" -> {10, 0},
"Velocity" -> {0, 1/20}|>,
<|"Mass" -> 1, "Position" -> {-10, 0}, "Velocity" -> {0, -1/20}|>},
400]


ParametricPlot[Evaluate[data[All, "Position", t]], {t, 0, 400}]

one gets

As a check, the separation at periapsis should be
$$r_p=a(1-e)=r_a\frac{1-e}{1+e}=20\frac{1-\frac{9}{10}}{1+\frac{9}{10}}=\frac{20}{19}\approx 1.05\tag{13}$$
and it seems to be.
A: You can use Mathematica to solve the numerical solution of this problem:
data  = NBodySimulation[
  "InverseSquare", {<|"Mass" -> 1, "Position" -> {0, 0}, 
    "Velocity" -> {0, .5}|>,
   <|"Mass" -> 1, "Position" -> {1, 1}, "Velocity" -> {0, -.5}|>}, 4]
ParametricPlot[Evaluate[data[All, "Position", t]], {t, 0, 4}]


