I want to create a program in $Mathematica$ that solves numerically the Three-body problem by Euler-Lagrange's equations. I was searching some methods to sucessfully do it. So I found a way to solve Two-body problem in http://www.maths.usyd.edu.au/u/joachimw/thesis.pdf (page 12). I found, also, how to minimize the number of generalized coordinates using relative-position vectors $\vec{s_i}=\vec{r_j}-\vec{r_k}$ and center of mass $\vec{r_G}$, as you can see it in the following scheme scanned from Goldstein, 3rd edition:
With this, I reached the Lagrangin to a Three-body system:
$$\mathcal{L}=\frac{1}{2}M\dot{r}_{G}^2+\frac{1}{2}\frac{m_1m_2}{M}|\vec{\dot{r}_1}-\vec{\dot{r}_2}|^2+\frac{1}{2}\frac{m_2m_3}{M}|\vec{\dot{r}_2}-\vec{\dot{r}_3}|^2+\frac{1}{2}\frac{m_1m_3}{M}|\vec{\dot{r}_1}-\vec{\dot{r}_3}|^2+2G \left[\frac{m_1m_2}{|\vec{{r}_1}-\vec{{r}_2}|}+\frac{m_2m_3}{|\vec{{r}_2}-\vec{{r}_3}|}+\frac{m_1m_3}{|\vec{{r}_1}-\vec{{r}_3}|}\right]$$
I want to use Euler-Lagrange's equations:
$$\frac{d}{dt}\left(\frac{\partial \mathcal{L}}{\partial{\dot{q}_i}}\right)-\frac{\partial \mathcal{L}}{\partial q_i}=0$$
where $q_i=r_G,s_1,s_2,s_3$.
Finding Lagrange's equations for $r_G$ was really easy, because $\frac{\partial \mathcal{L}}{\partial \dot{r}_{G}}$ is a conserved quantity ($r_G$ doesn't appear explicitly in $\mathcal{L}$). But finding Lagrange's equations for $s_i$, $i=1,2,3$ was a little bit confusing.
I have a doubt: is it true that $$\frac{d}{dt}\left(|\vec{{r}_j}-\vec{{r}_k}|\right)=|\vec{\dot{r}_j}-\vec{\dot{r}_k}|$$ or more specifically, $$\dot{s_i}=|\vec{\dot{r}_j}-\vec{\dot{r}_k}|~?$$ Because if it doesn't (and that's what I think) it wasn't useful using relative positions vectors. It would require new generalized coordinates as the $x,y,z$ components of $\vec{s_i}$.
What would be your approach to solve this problem? Could this be a good way to do it?