Mathieu eigenvalue problem in an interval

\[\begin{split} \mathbb{S} \begin{cases} \Omega = [0, L_x] \\ \textbf{k}=2\,\textbf{e}_x \\ u_{\text{D}}\vert_{\partial\Omega}=0 \\ u^{q=0}_{\text{e}, n}(x)=\sin(\sqrt{\lambda_{\text{e}, n}}x)~,~\lambda^{q=0}_{\text{e}, n}=\frac{n^2\pi^2}{L_x^2} \end{cases} \end{split}\]

import numpy as np
from dolfinx.fem import FunctionSpace

from lucifex.mesh import interval_mesh
from lucifex.fem import Constant
from lucifex.solver import evp, EigenvalueProblem, BoundaryConditions, OptionsSLEPc
from lucifex.viz import plot_line, create_cycler, save_figure
from lucifex.pde.eigen import mathieu


def create_solver(
    Lx: float,
    Nx: int,
    bcs: BoundaryConditions, 
    q: float,
    k: float,
    nev: int,
) -> EigenvalueProblem:
    mesh = interval_mesh(Lx, Nx)
    q = Constant(mesh, q)
    k = Constant(mesh, (k, ))

    fs = FunctionSpace(mesh, ('P', 1))
    slepc = OptionsSLEPc(
        eps_tol=1e-10,
        eps_target=0.0,
        eps_nev=nev,
        eps_ncv=50,
        eps_max_it=1000,
        eps_which='smallest_real',
    )
    u_solver = evp(mathieu, bcs, slepc)(fs, q, k)
    return u_solver


Lx = 1.0
Nx = 200
bcs = BoundaryConditions(
        ('dirichlet', lambda x: x[0], 0.0),
        ('dirichlet', lambda x: x[0] - Lx, 0.0),
    )
nev = 4
k = 2.0

q_axis = np.linspace(0, 20, num=10)
solvers = [create_solver(Lx, Nx, bcs, q, k, nev) for q in q_axis]

[s.solve() for s in solvers]
eigenvalues = [s.eigenvalues for s in solvers]
eigenreals = [[np.real(l) for l in lmbdas] for lmbdas in eigenvalues]
eigenreals_q0 = [(n * np.pi/ Lx) ** 2 for n in range(nev)]

lines = [
    (q_axis, [eigenreals[q_index][n] for q_index in range(len(q_axis))])
    for n in range(nev)
]
cyc = create_cycler(linestyle=['dashed'], marker=['o'], ms=[4], color=['black'])
slc = slice(1, None)
fig, ax = plot_line(
    lines[slc], 
    cyc=cyc,
    x_label='$q$', 
    y_label='$\lambda(q)$')
ax.set_xlim(q_axis[0], q_axis[-1])
ax.hlines(eigenreals_q0[slc], q_axis[0], q_axis[-1], color='black', lw=0.5, ls='-')
ax.set_yscale('log')
save_figure('eigenvalues', thumbnail=True)(fig)