Solution (source code)

= Solution

For a <normal mode> proportional to $e^{\mu t+ikx}$, put $\widehat D^2=D^2-k^2$. The three scalar equations become
$$
(\mu-\sigma\widehat D^2)\widehat\omega-\lambda D\widehat W=0,
$$
$$
(\mu-\sigma\widehat D^2)\widehat D^2\widehat W
+\lambda D\widehat\omega=-\sigma Ra\,k^2\widehat\theta,
$$
and
$$
(\mu-\widehat D^2)\widehat\theta=\widehat W.
$$
Multiplying through by the two scalar operators and eliminating $\widehat\omega,\widehat\theta$ yields
$$
\boxed{
\left[(\mu-\widehat D^2)(\mu-\sigma\widehat D^2)^2\widehat D^2
+\lambda^2D^2(\mu-\widehat D^2)
+\sigma Ra\,k^2(\mu-\sigma\widehat D^2)\right]\widehat W=0}.
$$