Suppose we have equation of the form:
$$H \Psi = E \Psi $$
where $H$ is Dirac Hamiltonian (also my question can be answered by people who are not familiar with Dirac Hamiltonian but familiar with numerical solution of different eigenvalue problems):
$$ H = \vec{\alpha} \vec{\nabla} + V(r) = H_0 + V(r)$$
and $\vec{\alpha}$ and $\beta$ are $4 \times 4$ Dirac matrices:
$$\vec{\alpha} = \left( \begin{matrix} 0 & \vec{\sigma} \\ -\vec{\sigma} & 0 \end{matrix} \right) \ , \ \ \beta = \left( \begin{matrix} 1 & 0 \\ 0 & -1 \end{matrix} \right)$$
$\vec{\sigma}$ are $2 \times 2$ Pauli matrices.
Here $H_0$ is free Dirac Hamiltonian (without potential, only kinetic term is preserved).
I think that the best way to start solving such an equation numerically is to introduce plane-wave basis (which diagonalizes free Hamiltonian $H_0$):
$$\Psi = \sum_n c_n \psi_n $$
Then we have a problem like:
$$H \sum_n c_n \psi_n = E \sum_n c_n \psi_n $$
My question is: what are the most commonly used ways to solve this equation numerically for arbitrary potential? Just to compute matrix elements $H_{nm}$ in this plane-wave basis, to find eigenvalue $E$ and coefficients $C_n$ from this computation? Or maybe there are some better methods?
