Solving this set of differential equations doesn't yield real eigenvalues I tried to solve a set of diff equations representing a 
$$M\ddot X=KX$$
$M=\begin{bmatrix}m & 0\\0 & m\end{bmatrix}$
$X=\begin{bmatrix}x_1' \\x_2'\end{bmatrix}$ (note that $x_1'=x_1-x_{1,eq}$ and $x_2'=x_2-x_{2,eq}$; this was done just to remove some constants from the equations)
$K=\begin{bmatrix}(8a^2b+k) & -k\\-k & (8a^2b+k)\end{bmatrix}$
I then assumed a solution of the form $x_1'(t) = A_1 e^{i\omega t}$ and $x_2'(t) = A_2 e^{i\omega t}$. I then plugged this into the set of differential equations and reduced to:
B$\begin{bmatrix}A_1\\A_2\end{bmatrix}=\begin{bmatrix}0\\0\end{bmatrix}$
where $B=\begin{bmatrix}(-m\omega^2-8a^2b-k) & k\\k & (-m\omega^2-8a^2b-k)\end{bmatrix}$
But when I do $det(B)=0$ (with $\omega$ as the eigenvalues), I get these as eigenvalues for $\omega$:
$$\omega_1 = -\frac{2ia\sqrt{2}\sqrt{b}}{\sqrt{m}}$$
$$\omega_2=\frac{2ia\sqrt{2}\sqrt{b}}{\sqrt{m}}$$
$$\omega_3=-\frac{\sqrt{2}\sqrt{-4a^2b-k}}{\sqrt{m}}$$
$$\omega_4=\frac{\sqrt{2}\sqrt{-4a^2b-k}}{\sqrt{m}}$$
These don't seem real-valued. These equations are supposed to represent a double well potential with a particle inside of each. Furthermore, the 2 particles are connected by a string. $k$ is the spring constant, $a,b$ are parameters of the potential function $=U(x)=b(x^2-a^2)^2$
 A: When I develop the system I get
$\begin{cases}
m\ddot x_1'=(8a^2b+k)x_1'-kx_2'=8a^2bx_1'+k(x_1'-x_2') \\
m\ddot x_2'=-kx_1'+(8a^2b+k)x_2'=8a^2bx_2'-k(x_1'-x_2') \\
\end{cases}\tag{E}$
So it is tempting to introduce two new variables $\begin{cases} u=x_1'+x_2' \\ v=x_1'-x_2' \end{cases}\iff\begin{cases} x_1'=\frac{u+v}{2} \\ x_2'=\frac{u-v}{2} \end{cases}$
By adding and substracting the lines of $(E)$ we now have : 
$\begin{cases} m\ddot u=8a^2bu\\ m\ddot v=8a^2bv+2kv \end{cases}\iff
\begin{cases} \ddot u-(\omega_1)^2u=0\\ \ddot v-(\omega_2)^2v=0 \end{cases}\quad\text{with}\quad\begin{cases} \omega_1=2a\sqrt{\frac{2b}{m}}\\ \omega_2=\sqrt{\omega_1^2+\frac{2k}{m}} \end{cases}$
I assumed all constants $a,b,m,k$ are positive from your physical description of the problem.
The solution in $(u,v)$ is :
$\begin{cases}
u=U_1\;e^{\omega_1 t}+U_2\;e^{-\omega_1 t} \\
v=V_1\;e^{\omega_2 t}+V_2\;e^{-\omega_2 t}
\end{cases}$
You have hyperbolic solutions, this is why when you assume $e^{i\omega t}$ form you get complex values for $\omega$. 
Plus $x_1,x_2$ are sums of these solutions, and you assumed only a single pulsation while you need two.
Depending of your initial conditions, it might be possible that the solutions reduce to $x_i=A\cosh(\omega_1t+\phi)\pm B\cosh(\omega_2t+\psi)$ (possibly with $\phi=\psi=0$). This is when $U1$ and $U_2$ have same sign, else it reduces to $\sinh$.
