Not to mention it’s diagonal plus convolution, so you can use iterative eigensolvers+FFTs+diagonal preconditioning to solve it at a cost linear in the number of unknowns. I imagine this comes from a band structure computation in quantum mechanics? You might be able to adapt GitHub - JuliaMolSim/DFTK.jl: Density-functional toolkit to your purposes (see eg Custom potential · DFTK.jl for how to incorporate a custom potential), although that might be overkill for this relatively simple problem.
2 Likes