Gauge transformation of differential equations I This is a follow-up question to Gauge transformation of differential equations. .
Let $y(x)$ be a solution to the following ODE:
\begin{eqnarray}
y^{''}(x) + a_1(x) y^{'}(x)+a_0(x) y(x)=0
\end{eqnarray}
Now define:
\begin{equation}
g(x):= \frac{y(x)+ r(x) y^{'}(x)}{r(x) \sqrt{a_0(x)} \exp(-1/2 \int a_1(x) dx)}
\end{equation}
where
\begin{equation}
r^{'}(x) + 1 - a_1(x) r(x)=0
\end{equation}
Then:
\begin{eqnarray}
&&g^{''}(x) + \\
&&\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \frac{1}{4} \left(\frac{2 a_0''(x)}{a_0(x)}+\frac{a_0'(x) \left(\frac{4}{r(x)}-2 a_1(x)\right)}{a_0(x)}-\frac{3 a_0'(x)^2}{a_0(x)^2}+4 a_0(x)+2
   a_1'(x)+\frac{8 a_1(x)}{r(x)}-a_1(x)^2-\frac{8}{r(x)^2}\right)g(x)=0
\end{eqnarray}
In[7]:= 
Clear[a0]; Clear[a1]; Clear[y]; Clear[r]; Clear[g]; Clear[m]; x =.; \
x0 =.;
r[x_] = Exp[Integrate[a1[x], x]] C[1] -  
   Exp[Integrate[a1[x], x]] Integrate[ Exp[-Integrate[a1[x], x]], x];
Simplify[r'[x] + 1 - a1[x] r[x]]
g[x_] = (y[x] + r[x] y'[x])/(
  r[x] Sqrt[a0[x]] Exp[-1/2 Integrate[a1[x], x]]);
Collect[(g''[x] + 
    1/4 (4 a0[x] + Derivative[1][a0][x]/a0[x] (4/r[x] - 2 a1[x]) - (
       3 Derivative[1][a0][x]^2)/a0[x]^2 + (
       2 (a0^\[Prime]\[Prime])[x])/a0[x] - a1[x]^2 + (8 a1[x])/r[x] + 
       2 Derivative[1][a1][x] - 8/r[x]^2) g[x]) //. {Derivative[2][y][
     x] :> -a1[x] y'[x] - a0[x] y[x], 
   Derivative[3][y][x] :> -a1'[x] y'[x] - a1[x] y''[x] - a0'[x] y[x] -
      a0[x] y'[x]}, {y[x], y'[x]}, Simplify]

Out[9]= 0

Out[11]= 0

Note that the result above can be used to generate ODEs whose solutions are known. For example let us take $j=1$ and  $B=C x_1$, $A=C x_1/x_2$ and :
\begin{eqnarray}
a_0(x)&=& (B C - A D)^2 \frac{x^{j-1}}{4(B+A x)^2 (B-D+(A-C) x)^2(D+C x)^2}\\
a_1(x)&=& \frac{2}{x}\\
\Longrightarrow\\
r(x)&=& \frac{x^2}{x_0} +x
\end{eqnarray}
then define:
\begin{eqnarray}
{\mathfrak P}_0&:=&x_0^2 x_2^2\\
{\mathfrak P}_1&:=&2 x_0 x_2 \left(x_2-4 C^2 x_1 (x_0 (x_1+x_2)-x_1 x_2)\right)\\
{\mathfrak P}_2&:=&x_2^2-8 C^2 x_0 \left(x_0
   \left(x_1^2+5 x_1 x_2+x_2^2\right)-x_1 x_2 (x_1+x_2)\right)\\
{\mathfrak P}_3&:=&-16 C^2 x_0 (2 x_0 (x_1+x_2)+x_1 x_2)\\
{\mathfrak P}_4&=&-8
   C^2 \left(3 x_0^2+3 x_0 (x_1+x_2)+x_1 x_2\right)\\
{\mathfrak P}_5&=&-8 C^2 (3 x_0+x_1+x_2)\\
{\mathfrak P}_6&=&-8 C^2
\end{eqnarray}
then we have:
\begin{equation}
g(x):= x\cdot \frac{y(x)+ r(x) y^{'}(x)}{r(x) \sqrt{a_0(x)}}
\end{equation}
Since from my answer to Looking for closed form solutions to linear ordinary differential equations with time dependent coefficients.  we know that $y(x)$ is expressed through hypergeometric functions we automaticaly know the solution to the following rather complicated ODE:
\begin{eqnarray}
g^{''}(x) + \left( \frac{\sum_{j=0}^6 {\mathfrak P}_j x^j}{4 C^2 x^2 (x+x_0)^2 (x+x_1)^2 (x+x_2)^2}\right) g(x)=0
\end{eqnarray}
Again my question in here would be find other cases where we can find close form solutions to ODEs which are too complicated to be handled using other methods.

You can think further about e.g. the effect the below approach apply on HEUN-type ODEs, or some superposition approaches with e.g. Solutions in terms of the hypergeometric functions etc.
Hopefully someone can get challenge on some quite advanced ODEs like:
$\dfrac{d^2u}{dr^2}+\left(\dfrac{1}{2(r+1)}+\dfrac{1}{2(r-1)}-\dfrac{1}{r^2}\right)\dfrac{du}{dr}-\dfrac{k_2}{2k_1^2}\left(\dfrac{1}{r+1}+\dfrac{1}{r-1}\right)u=0$
$\dfrac{d^2y}{ds^2}+\left(\dfrac{1}{2(s-6)}+\dfrac{1}{2(s+6)}-\dfrac{1}{s}\right)\dfrac{dy}{ds}+\left(\dfrac{6A-B}{2(s-6)}-\dfrac{6A+B}{2(s+6)}+A\right)y=0$ , $A\neq0$
$\dfrac{d^2f}{dr^2}+\left(\dfrac{2r}{r^2+1}-\dfrac{1}{\omega(r^2+1)^2}\right)\dfrac{df}{dr}-\dfrac{f}{\omega^2(r^2+1)^2}=0$
and so on
