Dirca's equation has is the following form:
$$i\hbar\gamma^{\mu}\partial_{\mu}\psi - mc\psi = 0$$
where $\mu = 0,1,2,3$ and $\gamma^{\mu}$ are the $4\times4$ matrices(in Dirac's representation):
\begin{align}
\gamma^0 =
\begin{bmatrix}
\mathbb{1} & 0 \\
0 & -\mathbb{1}
\end{bmatrix}
\gamma^k =
\begin{bmatrix}
0 & \sigma_k \\
-\sigma_k & 0
\end{bmatrix}
\end{align}
Where $\sigma_k$ are the Pauli Matrices ($k = 1,2,3$):
\begin{align}
\sigma_1 =
\begin{bmatrix}
0 & 1 \\
1 & 0
\end{bmatrix}
\sigma_2 =
\begin{bmatrix}
0 & -i \\
i & 0
\end{bmatrix}
\sigma_3 =
\begin{bmatrix}
1 & 0 \\
0 & -1
\end{bmatrix}
\end{align}.
Using the momentum e energy in quantum mechanics differential operators; $\vec{p} = -i\hbar\vec{\nabla}$, $E = ih\partial_t = i\hbar c\partial_0$ the Dirac's equation become:
$$(E\gamma^0 -c\gamma^{k} p_k - mc^2)\psi = 0$$
Breaking the four dimensional spinnor $\psi$ in two components the equation in a matrix for become:
\begin{align}
\begin{bmatrix}
(E - mc^2)\mathbb{1}& -\sigma_k p_k c\\
\sigma_k p_k c & -(E + mc^2)\mathbb{1}
\end{bmatrix}
\begin{bmatrix}
u\\
v
\end{bmatrix}
=0
\end{align}
To coupling the EM field the energy and the momentum change
\begin{align}
E \rightarrow E - e \Phi \\
\vec{p} \rightarrow \vec{p} - e \vec{A}
\end{align}
So the EM weak coupling Dirac's equation in matricial form become:
\begin{align}
\begin{bmatrix}
(E - e \Phi - mc^2)\mathbb{1}& -\sigma_k (p_k - e A_k) c\\
\sigma_k (p_k - e A_k) c & -(E - e \Phi + mc^2)\mathbb{1}
\end{bmatrix}
\begin{bmatrix}
u\\
v
\end{bmatrix}
=0
\end{align}
Now writing it in function of $u$ and $v$:
\begin{align}
(E - e \Phi - mc^2) u -(\sigma_k (p_k - e A_k) c)v =0 \\
\sigma_k (p_k - e A_k)c u -(E - e \Phi + mc^2)v =0
\end{align}
In the second equation the $u$ and $v$ relation become:
\begin{align}
\frac{\sigma_k (p_k - e A_k)c}{E - e \Phi + mc^2}u=v
\end{align}
For the non relativistic approach $e \Phi << mc^2$ and $E = mc^2$
So the $u$ and $v$ relation became:
\begin{align}
\frac{\sigma_k (p_k - e A_k)c}{2mc^2}u=v
\end{align}
Plugging this non relativistic relation in the first equation:
\begin{align}
(E - e \Phi - mc^2) u -\frac{(\sigma_k (p_k - e A_k)c)^2}{2mc^2}u =0
\end{align}
Renaming $E - mc^2= E_{NR}$ and reorganizing:
\begin{align}
\left( e \Phi +\frac{(\sigma_k (p_k - e A_k)c)^2}{2mc^2}\right)u =E_{NR}u
\end{align}
Focus in the $(\sigma_k (p_k - e A_k))^2$ term we have:
\begin{align}
(\sigma_k (p_k - e A_k))^2 &= (\sigma_i (p_i - e A_i)(\sigma_j (p_j - e A_j)\\
&= \sigma_i \sigma_j \Pi_i \Pi_j
\end{align}
Where $\Pi_i = p_i - e A_i$
From Pauli anti-commutation and commutation relations:
$$\sigma_i \sigma_j = \delta_{ij} + i \varepsilon_{ijk}\sigma_k $$
With that the quadratic term become:
\begin{align}
(\sigma_k (p_k - e A_k))^2 &= \Pi_i \Pi_i + i \varepsilon_{ijk}\sigma_k \Pi_i \Pi_j
\end{align}
The last term in $(\sigma_k (p_k - e A_k))^2 u$ is:
\begin{align}
i \varepsilon_{ijk}\sigma_k \Pi_i \Pi_j u &= i \varepsilon_{ijk}\sigma_k\left[(-i\hbar\partial_i - e A_i)(-i\hbar \partial_j - e A_j)\right] \\
&= i \varepsilon_{ijk}\sigma_k\left[-\hbar^2 \partial_i \partial_j + e^2 A_i A_j + i\hbar e(\partial_i A_j + A_i \partial_j)\right]u
\end{align}
The first two terms are symmetrical so they become:
\begin{align}
\varepsilon_{ijk} \partial_i \partial_j &= \frac{1}{2}\varepsilon_{ijk}(\partial_i \partial_j + \partial_j \partial_i) \\
&= \frac{1}{2}(\varepsilon_{ijk}\partial_i \partial_j +\varepsilon_{ijk}\partial_j \partial_i) \\
&= \frac{1}{2}(\varepsilon_{ijk}\partial_i \partial_j -\varepsilon_{jik}\partial_j \partial_i)\\
& = 0
\end{align}
The same thing are value for the $A_i A_j$ so the term become:
\begin{align}
i \varepsilon_{ijk}\sigma_k \Pi_i \Pi_j u &= -\hbar e \varepsilon_{ijk}\sigma_k\left[ \partial_i( A_j u) + A_i \partial_j (u)\right] \\
&= -\hbar e \varepsilon_{ijk}\sigma_k\left[ \partial_i( A_j)u+ A_j \partial_i (u) + A_i \partial_j (u)\right]
\end{align}
Where were applied the product rule in the first derivative. Now looking in the last two terms we can show that they vanish each other:
\begin{align}
\varepsilon_{ijk}\left[ A_j \partial_i + A_i \partial_j \right] &= \varepsilon_{ijk}A_j \partial_i + \varepsilon_{ijk}A_i \partial_j \\
&= \varepsilon_{ijk}A_j \partial_i - \varepsilon_{jik}A_i \partial_j \\
&= 0
\end{align}
The only term that does not vanish is the $k$-th component os the curl of $\vec{A}$
\begin{align}
i \varepsilon_{ijk}\sigma_k \Pi_i \Pi_j &= -\hbar e \varepsilon_{ijk}\sigma_k \partial_i( A_j) \\
&= -\hbar e \sigma_k \varepsilon_{kij}\partial_i A_j \\
&= -\hbar e \vec{\sigma}\cdot \vec{B}
\end{align}
With that we can write:
\begin{align}
(\sigma_k (p_k - e A_k))^2 &= \frac{(p_k - e A_k)^2}{2m} - \frac{\hbar e}{2m}\vec{\sigma}\cdot \vec{B}
\end{align}
Now finally putting all together the non relativistic Dirac's equation is:
\begin{align}
\left( e \Phi + \frac{(p_k - e A_k)^2}{2m} - \frac{\hbar e}{2m}\vec{\sigma}\cdot \vec{B}\right)u =E_{NR}u
\end{align}
