IMHO, the best approach to mass-spring systems is via matrix formulation,
$\underline{\underline{M}}\underline{\ddot{x}} + \underline{\underline{K}}\underline{x} = \underline{F} $,
where the vector $\underline{x}(t)$ collects the degrees of freedom of the problem, and finding the eigensolution of the homogenous system,
$\left[s^2 \underline{\underline{M}} + \underline{\underline{K}} \right]\underline{\hat{x}} = \underline{0} $.
In your example,
$\underline{\underline{M}} = \begin{bmatrix} m & 0 \\ 0 & m \end{bmatrix}$
$\underline{\underline{K}} = \begin{bmatrix} 2 k & -k \\ -k & 2 k \end{bmatrix}$
and the determinant of the matrix $s^2\underline{\underline{M}} +\underline{\underline{K}} $ reads
$\det (s^2\underline{\underline{M}} +\underline{\underline{K}} ) = (s^2 m + 2k)^2 - k^2 = m^2 s^4 + 4 km s^2 + 3 k^2$
the eigenvalues reads
$s_{1,2}^2 = - \left[ 2 \pm
1 \right]\dfrac{k}{m} $
and thus
$s_1^2 = - \dfrac{k}{m} $$\quad \rightarrow \quad$$s_{1\pm} = \pm j \sqrt{\dfrac{k}{m}}$
$s_2^2 = - 3 \dfrac{k}{m} $$\quad \rightarrow \quad$$s_{2\pm} = \pm j \sqrt{3\dfrac{k}{m}}$,
while the corresponding eigenvectors are:
$\underline{\hat{x}}_{1} =\begin{bmatrix} 1 \\ 1 \end{bmatrix}$,$\quad$$\underline{\hat{x}}_{2} =\begin{bmatrix} 1 \\ -1 \end{bmatrix}$.
The first eigensolution is the mode with the two masses oscillating "in phase" with the same amplitude and pulsation $\omega_1 = \sqrt{\frac{k}{m}}$.
The second eigensolution is the mode with the two masses oscillating "in anti-phase" with the same amplitude and pulsation $\omega_2 = \sqrt{3\frac{k}{m}}$.
We can use these eigensolutions to build the general solution of the homogeneous system as
$\underline{x}(t) = A_1 \underline{\hat{x}}_1 e^{j( \omega_1 t +\phi_1)} + A_2 \underline{\hat{x}}_2 e^{j( \omega_2 t+\phi_2)}$
$x_1(t) = A_1 \hat{x}_{11} e^{j( \omega_1 t +\phi_1)} + A_2 \hat{x}_{21} e^{j( \omega_2 t+\phi_2)}$
$x_2(t) = A_1 \hat{x}_{12} e^{j( \omega_1 t +\phi_1)} + A_2 \hat{x}_{22} e^{j( \omega_2 t+\phi_2)}$