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?

LogPlot[Abs@Rsol[0.5, 0.5][r], {r, 2.01, 50}]