I'm working on the solution of the equation
$$\tan^2u\partial^2_u y_2 + (2+\tan^2u)\tan u \partial_u y_2 -a^2\lambda_2y_2 - n^2(1+\cot^2u)y_2 = 0.$$ It is possible to write the above equation in terms of the Sturm-Liouville operator by just multiplying the equation by $\cos u$:
\begin{align*} \sin u \tan u \partial^2_u y_2 + (\cos u \tan u + \sec^2 u \sin u)\partial_u y_2 - \cos u(a^2\lambda_2 + n^2(1+\cot^2u))y_2 = 0\\ \implies \partial_u (\sin u \tan u \partial_u y_2) - \cos u (a^2\lambda_2 + n^2(1+\cot^2u))y_2 = 0. \end{align*} Now, let $\eta = \sin u$, hence,$$ \frac{\partial}{\partial u} = \frac{\partial\eta}{\partial u}\frac{\partial}{\partial\eta} = \cos u \frac{\partial}{\partial\eta}. $$ The last equation becomes \begin{align*} \cos u \partial_\eta \left(\eta^2\partial_\eta y_2\right) - \cos u (a^2\lambda_2y_2 + n^2(1+\cot^2u))y_2 &= 0\\ \implies \left(\eta^2\partial_\eta y_2\right) -\left(a^2\lambda_2 + n^2\left(1+\frac{\cos^2 u}{\sin^2 u}\right)\right)y_2 &= 0, \end{align*} but, since $\sin u = \eta,$ it implies that $\cos u = \sqrt{1-\eta^2}$ by trigonometric relations. Thus, \begin{align*} \partial_\eta \left(\eta^2\partial_\eta y_2\right) -\left(a^2\lambda_2 + n^2\left(1+\frac{1-\eta^2}{\eta^2}\right)\right)y_2 &= 0\\ \implies \partial_\eta \left(\eta^2\partial_\eta y_2\right) -\left(a^2\lambda_2 +\frac{n^2}{\eta^2}\right)y_2 &= 0 \end{align*}
Now, the problem is reduced to solve this last ODE, which I can see something similar to the associated Legendre ODE or Bessel ODE. The last one seems to be more plausible when we apply the change of variables $x=n/\eta$. By making this change, the equation becomes $$ x^2\partial^2_x y - (a^2\lambda_2 + x^2)y = 0$$
I want to reduce the equation to a known one to avoid the Frobenius method, but if it is not possible, I'm ok with that.
If possible, how can I write this last equation as a Bessel-type ODE?
Any help/idea of substitution will be very appreciated. Thanks in advance
You are almost done: the substitution $y(x)=\sqrt xz(x)$ in the equation \begin{equation} x^2\partial^2_x y(x) - (a^2\lambda_2 + x^2)y(x) = 0 \end{equation} gives \begin{equation} x^2\partial^2_x z(x)+x\partial_x z(x) - (a^2\lambda_2 +\frac14 +x^2)z(x)=0 \end{equation} which is a modified Bessel equation. The solutions are then a linear combination of the functions $K_{\sqrt{a^2\lambda_2+1/4}}(x)$ and $I_{\sqrt{a^2\lambda_2+1/4}}(x)$.
The idea of this substitution comes from the remark that the equation has exponential solutions for $x\to\infty$. If we expect that it is equivalent to a Bessel equation, the known asymptotic behavior of the Bessel functions ($x^{-1/2}\exp(\lambda x)$) suggests the change of function.