Asymptotic solution of a second-order ODE containing InverseFunction

Question

Essentially, I have a second-order differential equation given by ode below. In order to solve it, I need to obtain an asymptotic solution where $g(x)$ must vanish at infinity which will be used after as initial conditions (i.e. $g(\infty),g'(\infty)$) to evolve the differential equation using NDSolve. For the asymptotic solution, i have come up with the code below which uses the Frobenius method. The difficulty arises because the differential equation contains a function TP which is an InverseFunction of the rho[r,q]. I guess, my code does not work with function not user-defined. My questions are:

(1) Is there a way to fix my code? (2) Instead of using InverseFunction for inverting rho[r,q], is there a better way of inverting functions such that it works with the way the code is created?

asymp[p_] := {ode = (2/TP[x, q] - p/TP[x, q]^(p + 1)) TP[x, q]^2 dTP g'[x] + TP[x, q]^2 g''[x] - l (l + 1) g[x];
b = 1;
q = -1;
\[Alpha] = l + 1;
rho[r_, q_] := (2 b/(1 - q)) (1 - (b/r)^(1 - q))^(1/2) Hypergeometric2F1[1/2, 1 - 1/(q - 1), 3/2, 1 - (b/r)^(1 - q)]; 
TP = InverseFunction[Function[{r, q}, rho[r, q]], 1, 2];
dTP = (Sqrt[1 - (b/TP[x, q])^(1 - q)] ((b/TP[x, q])^(1 - q))^(1/(1 - q)) TP[x, q])/b;
g[r_] := Sum[a[i]/r^(i + \[Alpha]), {i, 0, ORDINF + 5}];
ORDINF = 3;
ss = FullSimplify[Series[ode, {x, \[Infinity], ORDINF}]];
eqsINF = Table[SeriesCoefficient[ss, i] == 0, {i, 2, ORDINF}];
yinf = Table[a[i], {i, 1, ORDINF - 1}];
seriesINF = Simplify[Solve[eqsINF, yinf]][[1]];
Rasymp = Collect[Simplify[Sum[a[i]/x^(i + \[Alpha]), {i, 0, ORDINF - 1}] /. seriesINF /. a[0] -> 1], x];
dRasymp = Collect[Simplify[D[Rasymp, x]], x]}

By the way, $p>0$ and $q<1$. Any help or hint would help. Thanks!

Answer 1

Here it is necessary to change the definition of the inverse function, and put the parameter q into the definition of the asymp[] along with the parameter l.

asymp[p_, q_, l_] := {
  b = 1;
  \[Alpha] = l + 1; ORDINF = 3;
  rho[r_] := (2 b/(1 - q)) (1 - (b/r)^(1 - q))^(1/
       2) Hypergeometric2F1[1/2, 1 - 1/(q - 1), 3/2, 
     1 - (b/r)^(1 - q)];
  TP[x_] := InverseFunction[rho][x];
  dTP[x_] := (Sqrt[
       1 - (b/TP[x])^(1 - q)] ((b/TP[x])^(1 - q))^(1/(1 - q)) TP[x])/
    b; ode = (2/TP[x] - p/TP[x]^(p + 1)) TP[x]^2 dTP [x] g'[x] + 
    TP[x]^2 g''[x] - l (l + 1) g[x];
  yinf = Table[a[i], {i, 1, ORDINF - 1}]; 
  g[r_] := Sum[a[i]/r^(i + \[Alpha]), {i, 0, ORDINF + 5}];

  ss = FullSimplify[Series[ode, {x, \[Infinity], ORDINF}]];
  eqsINF = Table[SeriesCoefficient[ss, i] == 0, {i, 2, ORDINF}];

  seriesINF = Simplify[Solve[eqsINF, yinf]][[1]];
  Rasymp = 
   Collect[Simplify[
     Sum[a[i]/x^(i + \[Alpha]), {i, 0, ORDINF - 1}] /. seriesINF /. 
      a[0] -> 1], x];
  dRasymp = Collect[Simplify[D[Rasymp, x]], x]}

Usage example

In[5]:= asymp[1, -1, 0] // Quiet

Out[5]= {1/(2 x^4) - 1/x^3 - 1/x^2}