Diagonalization of a matrix with operators as elements

How to diagonalize a hamiltonian matrix that has differential operators as elements? My matrix looks something like: $$\begin{bmatrix} A \frac{d^{2}}{{d\theta}^{2}}+ B_{1} & a\cos{(b\theta +c)}\\ a\cos{(b\theta +c)} & A \frac{d^{2}}{{d\theta}^{2}}+ B_{2} \end{bmatrix}$$ where A, $$B_{1}$$, $$B_{2}$$, a,b and c are all constants. How to do I find the spectrum of eigen values for such systems. Are there any standard methods available.

Let your matrix operate on a vector of functions (of $$\theta$$, in this case), and look for solutions where it produces a constant multiple of that vector. In other words, solve
$$$$\begin{bmatrix} A \dfrac{\mathrm d^{2}}{{\mathrm d\theta}^{2}} + B_{1} & & a\cos(b\theta +c) \\ a\cos{(b\theta +c)} && A \dfrac{\mathrm d^{2}}{{\mathrm d\theta}^{2}}+ B_{2} \end{bmatrix} \begin{bmatrix} f(\theta)\vphantom{\dfrac{\mathrm d^{2}}{{\mathrm d\theta}^{2}}} \\ g(\theta)\vphantom{\dfrac{\mathrm d^{2}}{{\mathrm d\theta}^{2}}} \end{bmatrix} =\lambda \begin{bmatrix} f(\theta)\vphantom{\dfrac{\mathrm d^{2}}{{\mathrm d\theta}^{2}}} \\ g(\theta)\vphantom{\dfrac{\mathrm d^{2}}{{\mathrm d\theta}^{2}}} \end{bmatrix} \tag{01}\label{eq01}$$$$ which is a coupled set of second-order differential equations for $$f(\theta)$$ and $$g(\theta)$$. To help solve this system and find values of $$\lambda$$ which produce periodic solutions so that there is no discontinuity in $$\theta$$, read about Mathieu functions.
• Sorry, I don't know anything about that transformation. But if you can find the operator's eigenvectors and eigenvalues, you can diagonalize it. Unfortunately, I think I'm wrong about the Mathieu functions... I originally didn't notice that the $B$'s had subscripts, and when I took $B_1=B_2=B$ and decoupled the equations by looking at $f+g$ and $f-g$ I got Mathieu functions. But for the more general $B_1\ne B_2$ I get fouth-order PDEs that I don't know how to solve analytically. Are you trying to get an analytic solution or a numerical solution? – G. Smith Nov 17 '18 at 5:28
• With Mathematica, I can numerically find the eigenvalues and eigenvector functions using the function NDEigensystem. For example, for $A=B2=B2=a=b=1, c=0$, I get that the first five eigenvalues are -0.2932, -0.2932, 0.405184, -1.28535, and -1.34277. – G. Smith Nov 17 '18 at 6:05