I'm trying to find the normal models of a particle with charge $q$ and mass $m$ in a $3$-dimensional harmonic oscillator potential with an applied uniform magnetic field $B=B_0 \hat{z}$. The potential has a natural frequency $\omega_0$. The Lagrangian is therefore
$$L = \frac{1}{2}m(\dot{x}^2 + \dot{y}^2 + \dot{z}^2) - \frac{1}{2}m\omega_0^2(x^2 + y^2 + z^2) - \frac{q B_0}{2c}(\dot{x}y - \dot{y}x)$$
The resulting equations of motion are \begin{align} &\ddot{x} = - \omega_0^2 x + \frac{q B_0}{mc} \dot{y} \\ &\ddot{y} = - \omega_0^2 y -\frac{q B_0}{mc} \dot{x} \\ &\ddot{z} = -\omega_0^2 z \end{align}
My understanding of finding the normal modes is that I should now fit these equations into the general form $m \ddot{\textbf{x}} = -k \textbf{x}$ for some matrix-valued $m$ and $k$; from there it's straightforward to find the normal modes as the eigenvectors of the matrix $-m^{-1}k$.
Unfortunately, the $x$ and $y$ equations of motion are not in a basic Hooke's law form -- they have a velocity term. I'm not sure how to proceed from the equations of motion, though I suspect I'm just missing something simple.