# How to resolve underflow occurrences?

I am solving a second order differential equation described by odey below. For the asymptotics, I have the following code which will be used as initial conditions for NDSolve.

(*asymptotics*)
asymp[p_,q_]:={
b0=1;
\[Alpha]=l+1;
ORDINF=5;
rho[r_]:=(2b0)/(1-q) Sqrt[1-(b0/r)^(1-q)]Hypergeometric2F1[1/2,1-1/(q-1),3/2,1-(b0/r)^(1-q)];
TP=InverseFunction[ConditionalExpression[rho[#1],#1>=0]&];
TPtab=Table[{r,TP[r]},{r,900,1000,1/10}];
fit=Normal[LinearModelFit[TPtab,y,y]];
g[r_]:=Sum[a[i]/r^(i+\[Alpha]),{i,0,ORDINF+5}];
ode= (2/fit-p/fit^(p+1))fit^2 D[fit,y] g'[y]+fit^2 g''[y]-l (l+1) g[y];
ss=FullSimplify[Series[ode,{y,\[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=Rationalize[Collect[Simplify[Sum[a[i]/y^(i+\[Alpha]){i,0,ORDINF-1}]/.seriesINF/.a[0]->1],y],rat],
dRasymp=Rationalize[Collect[Simplify[D[Rasymp,y]],y],rat]}


Notice that the asymptotic solution Rasymp and its derivative dRasymp have quite long expressions (with large numbers). Now I have the following initialization, and routine for solving the ODE. In the end, I have to obtain the object baremode[ell].

p = 0;
q = -1;
b0 = 1;
y0 = 1;
rules = {AccuracyGoal -> Infinity, PrecisionGoal -> 30, WorkingPrecision -> 30, MaxSteps -> 10000};
rat = 10^-30;
rho[r_] := (2 b0)/(1 - q) Sqrt[1 - (b0/r)^(1 - q)]Hypergeometric2F1[1/2, 1 - 1/(q - 1), 3/2, 1 - (b0/r)^(1 - q)];
ryy = InverseFunction[rho[#] &];
ry[y_] := Abs[ryy[y]]

(*differential equation*)
odey = (ry[y]^2/b0^2 R''[y] + (2/ry[y] - p/ry[y]^(p + 1)) ry[y]^2/b0^2 (Sqrt[1 - (b0/ry[y])^(1 - q)] Sign[y]) R'[y] - l (l + 1) R[y]) // Simplify(*nondimensionalized ODE*);

(*initial conditions for the scalar field (R[y], R'[y])*)
R1 = asymp[p, q][[1]];
R2 = asymp[p, q][[2]];
R[y0_, l0_] := R1 /. y -> y0 /. l -> l0
dR[y0_, l0_] := D[R2, y] /. y -> y0 /. l -> l0

yP = 10^3;
yM = -10^3;
Q = 1;
f[r_] := 1 - (b0/r)^(1 - q);
\[Psi][r_] := 1/r^p;
r0 = ry[y0];

For[ell = 0, ell <= 30, ell++, {R0p = Rationalize[R[yP, ell], rat];
dR0p = Rationalize[dR[yP, ell], rat];
R0m = Rationalize[R[yM, ell], rat];
dR0m = Rationalize[dR[yM, ell], rat];

BCp = {R[yP] == R0p, R'[yP] == dR0p};
EQp = {(odey /. l -> ell) == 0};
Rsolp = NDSolveValue[Union[EQp, BCp], R, {y, yP, y0}, rules,
Method -> "StiffnessSwitching"];

BCm = {R[yM] == R0m, R'[yM] == dR0m};
EQm = {(odey /. l -> ell) == 0};
Rsolm = NDSolveValue[Union[EQm, BCm], R, {y, yM, y0}, rules,
Method -> "StiffnessSwitching"];

rp = Rationalize[ Rsolp[y0], rat];
rm = Rationalize[ Rsolm[y0], rat];
drp = Rationalize[ Rsolp'[y0], rat];
drm = Rationalize[ Rsolm'[y0], rat];
s = Rationalize[-Q b0 Sqrt[4 Pi (2 ell + 1)]/ry[y0]^2, rat];
c1f = (rm s)/(drp rm - drm rp);
r0 = ry[y0];
baremode[ell] = Sqrt[(2 ell + 1)/(4 Pi)] 1/b0 (b0/f[r0]  Sqrt[1 - (b0/ry[Abs[y0]])^(1 - q)]) c1f Rsolp'[y0];}]


But if I run the code, I get an error message that says an Underflow occurred in computation. I understand that along the numerical calculation, it encounters an extremely small numerical value that cannot be handled with the precision number set by MMA. I tried to track down where does the small number appear in my calculation, but I have not seen it yet. Moreover, I was having a hard time dealing with InverseFunction. I suspect, the problem is partly caused by my bad implementation. I also tried increasing the WorkingPrecision but still the problem still appears. Maybe someone could help.