Two-mass spring system in $x$-$y$ plane motion I wrote a coupled differential system that describes two masses coupled by one spring.  Both masses are free to move in x-y plane.  To test the model, I plotted a position of two masses in x-coordinate under three initial conditions:

*

*for mass 1 $x_0=A_o$ and $y_0=0$, for mass 2 $x_0= +A_o$ and $y_0=0$

*for mass 1 $x_0=A_o$ and $y_0=0$, for mass 2 $x_0= -A_o$ and $y_0=0$

*for mass 1 $x_0=A_o$ and $y_0=0$, for mass 2 $x_0= 0 $ and $y_0=0$
where $\omega=1$ and $k=0.1$
The plots for each should be as follows:

*

*two cosine waves in phase

*two cosine waves out of phase by 90 degree or cosine and sine

*two FM modulated waves out of phase ( two waves with beats)

example of one wave with beats

The visual model of the system is shown below,

My equations are
$  \frac{d^2x_1}{dt^2}+\omega_1^2x_1+k(x_1-x_2)\frac{x_1}{\sqrt{x_1^2+y_1^2}}=0   $
$  \frac{d^2y_1}{dt^2}+\omega_1^2y_1+k(y_1-y_2)\frac{y_1}{\sqrt{x_1^2+y_1^2}}=0   $
$  \frac{d^2x_2}{dt^2}+\omega_2^2x_2+k(x_2-x_1)\frac{x_2}{\sqrt{x_2^2+y_2^2}}=0   $
$  \frac{d^2y_2}{dt^2}+\omega_2^2y_2+k(y_2-y_1)\frac{y_2}{\sqrt{x_2^2+y_2^2}}=0   $
When I test the above equation, I get the correct plot for condition 1 and 2.  But for condition 3 I get division by zero error.  To address the issue, I used 0.000000001 instead of 0.  Then the  plot is inconsistent, there are not beats.  If I had the model made only in x-coordinate only, there would be no last term (multiplier with numerator and denominator). The system then gives all 3 conditions plots correctly.  Looks like my model is not correct.  Whats wrong with my model?
 A: I will change the coordinate to the absolute coordinate related to the origin:
\begin{array}
\text{ Define }\, \text{notions:}&\\
   \text{position of  mass 1 } & \vec{r}_1 =  \left( x_1(t) , y_1(t)\right) ;\\
   \text{position of mass 2 } & \vec{r}_2 = \left( x_2(t) , y_2(t)\right);\\
\text{ equilibrium position} &  \vec{r}_{0} =\left( x_{20}-x_{10} , y_{20}-y_{10}\right);\\
\text{ equilibrium length } &  {\ell}_{0} =\sqrt{ \left(x_{20}-x_{10}\right)^2 + \left(y_{20}-y_{10}\right)^2};\\
 \text{relative position }   &   \vec{r}_{12} =\left( x_2(t)-x_1(t) , y_2(t)-y_1{t}\right); \\
\text{ length } &  {\ell}(t) =\sqrt{ \left(x_{2}(t)-x_{1}(t)\right)^2 + \left(y_{2}(t)-y_{1}(t)\right)^2};\\
 \text{Force of mass 1 }   &   \vec{F}(t) = k \frac{\left( x_2(t)-x_1(t) ,\, y_2(t)-y_1(t)\right)}{\ell)t)} \left( \ell(t) - \ell_0\right) = k \left( x_2(t)-x_1(t) ,\, y_2(t)-y_1{t}\right) \left( 1 - \frac{\ell_0}{\ell(t)}\right) ; \\
 \text{Force component x }   &   F_x(t) = k \left( x_2(t)-x_1(t)\right) \left( 1 - \frac{\ell_0}{\ell(t)}\right); \\
 \text{Force component y }   &   F_x(t) = k \left( y_2(t)-y_1(t)\right) \left( 1 - \frac{\ell_0}{\ell(t)}\right); \\
\end{array}
