Solution (source code)

= Solution

Apply $i\gamma^\nu\partial_\nu+m$ to the <Dirac equation>. Since partial derivatives commute, only the symmetric part of the <gamma matrix> product contributes:
$$
\begin{aligned}
0&=(i\gamma^\nu\partial_\nu+m)(i\gamma^\mu\partial_\mu-m)\psi\\
&=-\left(\tfrac12\{\gamma^\nu,\gamma^\mu\}\partial_\nu\partial_\mu+m^2\right)\psi
=-(\Box+m^2)\psi.
\end{aligned}
$$
Here $\Box=\partial_t^2-\nabla^2$ in the signature fixed above. The remaining operator is scalar in spinor indices, giving
$$
\boxed{(\Box+m^2)\psi_\alpha=0\quad\text{for each of the four components}.}
$$
This is the <Klein-Gordon equation> for each component, although the first-order <Dirac equation> additionally couples them.