I first thought about solving the eigenvalues and position representation eigenfunctions by assuming the eigenfunctions are of product form ##\psi(x,y) = \alpha(x)\beta(y)## and using separation of variables to get ODE:s for ##\alpha## and ##\beta##, but if the coefficients ##\sigma_{i}## are matrices instead of numbers, that probably makes the situation more complicated...