# Numerical solution of an ODE using NDSolve

I am to find the solution $R(\rho)$ to the following differential equation,

$\frac{\partial}{\partial\rho}\left(A^2\frac{\partial}{\partial\rho}R(\rho)\right)-l(l+1)R(\rho)=0$

where $A$ is a numerically inverted function of $\rho$. My code is,

 rho[r_, b_, 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, b, q}, rho[r, b, q]], 1, 3];
ode = {D[(Tp[\[Rho], 0, -1])^2 D[R[\[Rho]], \[Rho]], \[Rho]] - l (l + 1)R[ [Rho]] + 2 \[Xi] a^2/(\[Rho]^2 + a^2) R[\[Rho]] == 0};
Bc = {R[0] == 1, R'[0] == 1};
rules = {AccuracyGoal -> Infinity, PrecisionGoal -> 20, WorkingPrecision -> 30, MaxSteps -> 10000};
Rsol = R /. First@NDSolve[Union[ode, Bc], {R}, {\[Rho], 0, 10}, rules]


This is what Mathematica says: I don't know what is wrong with my code.

It is necessary to express the derivative of the inverse function in terms of the derivative of the original function using the rule $dy/dx=1/dx/dy$:

f = InverseFunction[
Function[{y, b,
q}, (2 b/(1 - q)) (1 - (b/y)^(1 - q))^(1/2) Hypergeometric2F1[
1/2, 1 - 1/(q - 1), 3/2, 1 - (b/y)^(1 - q)]], 1, 3];

f1[x_, b_, q_] := (
Sqrt[1 - (b/f[x, b, q])^(1 - q)] ((b/f[x, b, q])^(1 - q))^(1/(1 - q))
f[x, b, q])/b
ode = 2 f[x, b, q]*f1[x, b, q]*g'[x] + f[x, b, q]^2*g''[x] -
l (l + 1) g[x] == 0;
b = 1; q = -1; sol =
ParametricNDSolveValue[{ode, g[1] == 1, g'[1] == 1},
g, {x, 1, 10}, {l}]
Plot[Table[sol[l][x], {l, 0, 5, 1}], {x, 1, 10}]


Note that in a certain range of parameters the solution becomes complex, so you need to print the real part or the absolute value of the function

• Hi @Alex. I cannot reproduce your expression for f1[x,b,q]. How did you arrive with that? Commented Sep 3, 2018 at 5:05
• I'm not quite sure how did you come up with the expression for f1. But I used 1/(dx/dy) and I obtained different expression. Commented Sep 3, 2018 at 7:35
• `rho[x_, b_, q_] := (2 b/(1 - q)) (1 - (b/x)^(1 - q))^(1/2) Hypergeometric2F1[ 1/2, 1 - 1/(q - 1), 3/2, 1 - (b/x)^(1 - q)]; 1/D[rho[x, b, q], x] // FullSimplify (Sqrt[1 - (b/x)^(1 - q)] ((b/x)^(1 - q))^(1/(1 - q)) x)/b' Commented Sep 3, 2018 at 7:57
• Right. There's a problem with the code. I can't successively run different values of $q$. For example, I run for $q=-1$ and it spits out a plot. Then, I change the value of $q$ to $-2$, and error messages appear. What can I do? Commented Sep 3, 2018 at 8:13
• I checked, that with q = -2 it works. Each time you need to restart the Mathematica for each new task. Make a list of parameters, we will try to rewrite the code and run the task in one session. Commented Sep 3, 2018 at 8:30