Simple method for calculating and plotting the band structure of a known periodic potential?

A thousand pardons if this is trivial, but I’ve been stuck here for hours.

I’m trying to compute the spectrum and eigenfunctions (i.e. band structure) of the eigenvalue equation $u” + k^2epsilonleft(xright)u = k^2q^2u$, where k is a known parameter, $epsilonleft(xright)$ is a known periodic function and $q^2$ is essentially the eigenvalue that I need to solve for and plot.

Now, I know about Bloch’s theorem, and substituting the Bloch ansatz gives another ODE whose solution has to be periodic, but couldn’t figure out how to proceed from there.

Most approximation methods that I found are in the context of condensed matter physics and seemed fairly involved, but I am working in the context of photonic band structures. Literature on photonic band structures of this sort don’t present any intermediate computation steps, like it’s something very standard.

How can I go about solving this? Also, are there ODE solvers that can handle spectrum calculations for periodic coefficients or with constraints of periodic solutions?