# Asymptotic solution of a second-order ODE containing InverseFunction

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!

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}