# Statement of the problem

I need to (numerically) solve an eigenproblem of the type

$$-\omega^2\mathcal{D}_1\vec{x}=\mathcal{D}_2\vec{x}$$ on the interval $[-1,1]$, where $\mathcal{D}_1$ and $\mathcal{D}_2$ are differential operators with a generalised boundary condition of the form

$$\omega^2\mathcal{D}_3\vec{x}(±1)+a~\omega\mathcal{D}_4\vec{x}(±1)+b~\mathcal{D}_5\vec{x}(±1)=0~.$$

where $\mathcal{D}_3$, $\mathcal{D}_4$ and $\mathcal{D}_5$ are more differential operators.

The critical feature here being the fact that the eigenvalue appears in the boundary condition.

# My personal take on it

SLEPc offers a wide range of method for dealing with generalised eigenvalue problems including the so-called polynomial eigenvalue problem. I used a spectral representation to cast my differential problem into a regular quadratic matrix eigenvalue problem of the form

$$(A\omega^2+B\omega+C)\vec{x}=0$$

I have used row replacement to include the boundary condition into my matrices. Since the first power of $\omega$ only appears in the boundary condition, the $B$ matrix is mostly empty.

For clarity, here is an idea of what the matrices look like

I expect that, if everything is done correctly, feeding this system into a dedicated solver for quadratic eigenvalue problems should spit out the expected eigenvalues (in some easy case, my particular system has analytical solutions). However, I haven't been able to recover these.

# My question

I have already spent a decent amount of time trying to find bugs in my code and will continue to do so. But in the meantime, I figured I might ask if anyone had a similar problem and had a different approach. I am especially anxious to know if perhaps there are some results unknown to me that might render my whole approach too naive and hopeless.

Any help would be greatly appreciated.