We consider a simple cubic crystal with one atom in the basis. The primitive lattice vectors are,

\begin{equation} \vec{a}_1=a\,\hat{x},\quad \vec{a}_2=a\,\hat{y},\quad\vec{a}_3=a\,\hat{z}. \end{equation}

Each atom has 6 nearest neighbors and 12 next-nearest neighbors. Let $u^x_{l,m,n}$ the displacement of an atom in the unit cell at position $\vec{r}=l\vec{a}_1+m\vec{a}_2+n\vec{a}_3$ in the $x$-direction from its equilibrium position. If the restoring force that pushes the atoms back to their equilibrium positions is modelled with linear springs with spring constants $C_{1}$ for nearest neighbors and $C_{2}$ for next-nearest neighbors, then the equations of motion for the atoms are,

\begin{equation} \begin{split} M \frac {d^{2}u^{x}_{lmn}}{dt^{2}} & = C_1\left( \left(u^{x}_{(l+1)mn} - u^{x}_{lmn}\right)+\left(u^{x}_{(l-1)mn} - u^{x}_{lmn}\right)\right)\\ &+\frac{C_2}{2} \Big[\left(u^{x}_{(l+1)(m+1)n} - u^{x}_{lmn}\right)+\left(u^{x}_{(l+1)(m-1)n} - u^{x}_{lmn}\right)+ \left(u^{x}_{(l-1)(m+1)n} - u^{x}_{lmn}\right)+ \left(u^{x}_{(l-1)(m-1)n} - u^{x}_{lmn}\right)\\ &+\left(u^{x}_{(l+1)m(n+1)} - u^{x}_{lmn}\right)+\left(u^{x}_{(l+1)m(n-1)} - u^{x}_{lmn}\right)+ \left(u^{x}_{(l-1)m(n+1)} - u^{x}_{lmn}\right)+ \left(u^{x}_{(l-1)m(n-1)} - u^{x}_{lmn}\right)\\ &+\left(u^{y}_{(l+1)(m+1)n} - u^{y}_{lmn}\right)+\left(u^{y}_{(l+1)(m-1)n} - u^{y}_{lmn}\right)+ \left(u^{y}_{(l-1)(m+1)n} - u^{y}_{lmn}\right)+ \left(u^{y}_{(l-1)(m-1)n} - u^{y}_{lmn}\right)\\ &+\left(u^{z}_{(l+1)m(n+1)} - u^{z}_{lmn}\right)+\left(u^{z}_{(l+1)m(n-1)} - u^{z}_{lmn}\right)+ \left(u^{z}_{(l-1)m(n+1)} - u^{z}_{lmn}\right)+ \left(u^{z}_{(l-1)m(n-1)} - u^{z}_{lmn}\right)\Big] \end{split} \end{equation}

There are similar equations for the displacements in the $y-$ and $z-$directions. The normal modes must be eigenfunctions of the symmetry of the system. Because of the translational symmetry of the crystal, the normal modes must have the form,

\begin{equation} \vec{u}_{l,m,n}=\vec{u}_{\vec{k}}\exp\left( i(lk_xa+mk_ya+nk_za)\right). \end{equation}

Substituting this form for the normal modes into the equations of motion above results in a set of algebraic equations that can be written in matrix form,

$$\frac{M}{C_1}\omega^2\vec{u}_{\vec{k}}= \textbf{m} \vec{u}_{\vec{k}},$$

where the elements of the matrix $\textbf{m}$ are:

$$m_{11} = 2\left(1-\cos(k_xa)\right)+2\frac{C_2}{C_1}\left(2-\cos(k_xa)\cos(k_ya)-\cos(k_xa)\cos(k_za)\right)$$ $$m_{12} = m_{21} = 2\frac{C_2}{C_1}\sin(k_xa)\sin(k_ya)$$ $$m_{13} = m_{31} = 2\frac{C_2}{C_1}\sin(k_xa)\sin(k_za)$$ $$m_{22} = 2\left(1-\cos(k_ya)\right)+2\frac{C_2}{C_1}\left(2-\cos(k_xa)\cos(k_ya)-\cos(k_ya)\cos(k_za)\right)$$ $$m_{23} = m_{32} = 2\frac{C_2}{C_1}\sin(k_ya)\sin(k_za)$$ $$m_{33} = 2\left(1-\cos(k_za)\right)+2\frac{C_2}{C_1}\left(2-\cos(k_xa)\cos(k_za)-\cos(k_ya)\cos(k_za)\right)$$

The eigenvalues $\lambda = \frac{M\omega^2}{C_1}$ can be found by setting the determinant equal to zero, $|\textbf{m} -\lambda \textbf{I}| = 0$.

This can be written as a cubic equation in $\lambda$.

$$-\lambda^3 +(m_{11}+m_{22}+m_{33})\lambda^2+(m_{12}^2+m_{13}^2+m_{23}^2 - m_{11}m_{22} - m_{11}m_{33} - m_{22}m_{33})\lambda + m_{11}m_{22}m_{33} +2m_{12}m_{13}m_{23} - m_{12}^2m_{33} - m_{13}^2m_{22}-m_{23}^2m_{11} =0$$

Cubic equations can be solved using Cardano's formula. The standard form of a cubic equation is,

$$a\lambda^3 + b\lambda^2 + c\lambda + d =0$$

Code that will calculate the coefficients $a,b,c,d$ and the matrix elements from the $\vec{k}$-vector is,

$k_xa=$  $k_ya=$  $k_za=$   $\frac{C_2}{C_1} = $  

$\sqrt{\frac{M}{C_1}}\omega = $

$\Gamma: 0,0,0$

$X: \pi/a,0,0$

$M: \pi/a,\pi/a,0$

$R: \pi/a,\pi/a,\pi/a$