I am playing with the one dimensional massless Dirac equation $$-i\partial_x \sigma_x \Psi(x) = \epsilon \Psi(x) \, .$$ Solving it analytically with some boundary conditions is an easy task. I am more interested in the numerical approach. One can rewrite the equation in matrix form using finite symmetric difference and obtain.
$$M\Psi = V$$ $$M=\begin{pmatrix} 1 & 0 & 0 & \cdots & 0 & 0 \\ -t & -\epsilon & t & \cdots & 0 & 0 \\ 0 & -t & -\epsilon & \ddots & 0 & 0 \\ \vdots & \vdots & \ddots & \ddots & t & 0 \\ 0 & 0 & 0 & -t & -\epsilon & t \\ 0 & 0 & 0 & 0 & 0 & 1 \\ \end{pmatrix}$$ $$t = -i\sigma_x/2a$$ $$\Psi=(\Psi_1, \Psi_2, \dots, \Psi_N)^{T}$$ $$V=(\Psi_0, 0, \dots, 0, \Psi_L)^{T}$$ Where the boundary conditions are implemented $\Psi(0) = \Psi_0$ and $\Psi(L)=\Psi_L$. Here I am stuck. How can I proceed to calculate the energies and corresponding wavefunction?