NDSolve solution plugged back into ODE quickly diverges

I have the following ODE with complex $$\omega$$:

$$\qquad f(r)\frac{d}{dr}(f(r)\frac{dR}{dr})+(\omega^2-f(r)(\frac{1}{r^2}+\mu^2))R(r)=0,$$

and with $$f(r)=1-2/r$$. I want to solve it starting at $$r=2$$ outwards for a specific set of $$\{\omega,\mu\}$$. For that I use NDSolve and the code

f[r_] := 1 - 2/r (*M=1*)
Diff = ((f[r] D[f[r] D[#, r],r] + (ω^2 - f[r] (1/r^2+μ^2)) #)) &;
eq = ((Diff @ R[r]) //. {ω -> ωnr + I ωni, μ->0.01} // Simplify) == 0;
sol =
ParametricNDSolve[
{eq, R[2.01] == 1 + I, R'[2.01] == 1 + I},
R, {r, 2.01, 100}, {ωnr, ωni},
MaxStepSize -> 0.001];
Rsol = R /. sol;
Shouldbezero =
Diff@(Rsol[0.5, 0.5][r]) //. {ω -> 0.5 + I 0.5, μ ->0.01} /. {M -> 1} // Simplify;
LogPlot[Shouldbezero // Abs, {r, 2.01, 50}]

As can be seen above, I plugged the solution that NDSolve provides back into the ODE, to check, whether the result is actually a solution. To my surprise, the numerical solution seems to be off, by quite a bit.

Plot of the absolute value of R(r) plugged back into the ODE: I have played around with AccuracyGoal, etc., but could only achieve minimal improvement. How come Mathematica gives me such a poor result? Or is this due to the nature of the equation?

• I don't think it's a poor result. Do notice the order of magnitude of solution is 14 around r==50. – xzczd Mar 22 at 4:19
• Looks like an issue of relative vs. absolute error. – J. M. will be back soon Mar 22 at 12:06 