# How to solve this 2nd-order ODE with quadratic coefficients?

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 solvable by Mathematica?

I tried the following to reduce it to a 2nd-order ODE of $$b$$ with coefficients quadratic in $$y$$. But DSolve only gives a useless DifferentialRoot form after a long wait.

variables = {a[y], b[y]};
Fop1[F_, pm_] :=
(r y + s) F - I (p D[F, y] + q (y D[F, y])) +
pm (-I) ((w y + t) F - I (u D[F, y] + v (y D[F, y])));
lhs = {Fop1[b[y], 1], Fop1[a[y], -1]};
eqe =
FullSimplify[
Eliminate[
Flatten[{D[(lhs - λ variables) // First, y],
lhs - λ variables}] == 0],
{a[y], a'[y]}]]
DSolve[eqe, b[y], y]


However, when $$u,\,v=0$$ or $$p,\,q=0$$, it becomes solvable, although the coefficients are still quadratic polynomials of $$y$$. Therefore, I was wondering if the more general case could be tackled as well. But I don't know how to proceed.

• Looks very unlikely to get a general solution to the underlying ODE with quadratic coefficients. If you play around with general equations of the form it requires several coefficients to be zero to get a special function representation to be found. – KraZug Jan 21 '19 at 10:47

Think that

$$(d_1-i d_2)b(y) = \lambda a(y)\\ (d_1+i d_2)a(y) = \lambda b(y)$$

can be handled as

$$(d_1+i d_2)(d_1-i d_2)b(y) = \lambda (d_1+i d_2) a(y) = \lambda^2 b(y)$$

then

$$(d_1^2+i(d_2d_1-d_1d_2)+d_2^2)b(y) = \lambda^2 b(y)$$

In this case $$d_1d_2 \ne d_2d_1$$ and

$$d_1^2 +i(d_2d1-d_1d_2)+ d_2^2 = p_2(y)\partial_x^2+p_1(y)\partial_x+p_0(y)$$

with

$$\left\{ \begin{array}{rcl} p_2(y) & = &-(p+q y)^2-(u+v y)^2\\ p_1(y) & = &-p (q+i (2 r y+2 s-v))-q^2 y-i q (2 y (r y+s)+u)-(u+v y) (2 i t+v+2 i w y)\\ p_0(y) & = & -i r (p+q y)-w (p+q y)+(r y+s)^2+r (u+v y)+(t+w y)^2-i w (u+v y) \end{array} \right.$$

and finally

$$p_2(y)b''(y)+p_1(y)b'(y)+(p_0(y)-\lambda^2)b(y) = 0$$

This DE can be solved with a series expansion solution, proposing

$$b(y) = \sum_{k=0}^n\alpha_k y^k$$

Clear[p, q, u, v, r, s, w, t]
d1[f_, y_] := -I (p + q y) D[f, y] + (r y + s) f
d2[f_, y_] := -I (u + v y) D[f, y] + (w y + t) f
d11 = d1[d1[b[y], y], y]
d22 = d2[d2[b[y], y], y]
d12 = d2[d1[b[y], y], y] - d1[d2[b[y], y], y]
DE = d1[d1[b[y], y], y] + I d1d2 + d2[d2[b[y], y], y] - lambda^2 b[y]
c0 = D[DE, b[y]]
c1 = D[DE, b'[y]]
c2 = D[DE, b''[y]]
Operator = c2 D[#, {y, 2}] + c1 D[#, y] + c0 # &;
n = 4;
Sumb = Sum[Subscript[alpha, k] y^k, {k, 0, n}];
res = Operator[Sumb] /. {Subscript[alpha, 0] -> Subscript[b, 0],
Subscript[alpha, 1] -> Subscript[b, 1]};
coefs = CoefficientList[res, y];

For[k = 1; Alphas = {}; equsk = equs[[1]]; subs = {},
k <= Length[equs] - 1, k++,
solalphak = Solve[equsk, Subscript[alpha, k + 1]];
AppendTo[Alphas, Subscript[alpha, k + 1] /. solalphak][[1]];
AppendTo[subs, solalphak];
equsk = equs[[k + 1]] /. Flatten[subs]
]
Alphas


So we can extract the $$\alpha_k$$ in Alphas

• Thanks. But it's wrong from the 4th line as you can easily check using my code. – xiaohuamao Jan 21 '19 at 17:16
• @xiaohuamao Yes. A forgotten # in the Operator definition. Now it is fixed. Thanks. – Cesareo Jan 21 '19 at 18:12
• I mean your Latex derivation is wrong from the 4th line. – xiaohuamao Jan 21 '19 at 18:22
• @xiaohuamao Yes. A wrong defined operator. Now I hope it is fixed. Thanks. – Cesareo Jan 21 '19 at 19:31
• The expansion is an undetermined Laurent series with unknown negative powers. Not as simple as you've assumed. And the 1st-order operators are non-Abelian. And I still find sign errors in your 2nd-order operator. – xiaohuamao Jan 21 '19 at 19:48