Consider an ODE eigensystem $$\begin{bmatrix} 0 & d_1-\mathrm id_2 \\ d_1+\mathrm id_2 & 0 \end{bmatrix} \begin{bmatrix} a(y) \\ b(y) \end{bmatrix} = \lambda \begin{bmatrix} a(y) \\ b(y) \end{bmatrix}, $$ where $$d_1=-\mathrm i(p+qy)\partial_y+ry+s$$ $$d_2=-\mathrm i(u+vy)\partial_y+wy+t,$$ $p,q,r,s,u,v,w,t$ are just real constants, and $\mathrm i$ is the imaginary unit. Is it analytically solvable?
I reduce it to a 2nd-order ODE of $b$ with coefficients quadratic in $y$ $$\alpha b''(y) + \beta b'(y) + \gamma b(y)=-\lambda^2 b(y)$$ where $$\alpha=(p+q y)^2+(u+v y)^2$$ $$\beta=p (q+2 i s-i v)+u (v+iq+2 it)+(2 i p r+q^2+2 i q s+2 i t v+2 i u w+v^2)y+2 i (q r+v w)y^2 $$ $$\gamma=-s^2-t^2+(p+i u) (w+i r)+[w (q-2 t+i v)-r (-i q+2 s+v)]y-(r^2+w^2)y^2 $$ When $u,v=0$ or $p,q=0$, it is solvable, although the coefficients are still quadratic polynomials of $y$. I was wondering if the more general case could be tackled as well. But I don't know how to proceed.