Finding eigenvalues of a differential operator

I am trying to get the eigenvalues of the following differential operator

$$L\psi(r) = -f \partial_r (f \partial_r \psi(r)) + V \psi(r)$$

which must satisfy (obviously)

$$L \psi(r) = \omega^2 \psi(r)$$

where I want to acquire both the real and imaginary part of $\omega$. To my problem, we have

$$f = 1 - \frac{2M}{r}$$ and $$V = f \left( \frac{l (l-1)}{r^2} + \frac{2 (1-S^2) M}{r^3} \right)$$ with $M = 1, l = 2, S=2$. The boundary conditions are $\psi(2M) = 0, \psi(\inf) = 0$. (To whomever may care, I am getting the real and imaginary oscillations of a Schwarzschild black hole. It is well studied in the literature, but I need to recover the result).

I tried three different ways to do it:

1) Using NDEigensystem:

f = 1 - 2*(M/r); 
V = f*(l*((l - 1)/r^2) + (2*(1 - S^2)*M)/r^3); 

M = 1; l = 2; S = 2; 
\[ScriptCapitalL] = f*D[f*D[\[Chi][r], r], r] + V*\[Chi][r]; 
\[ScriptCapitalB] = DirichletCo - ndition[\[Chi][r] == 0, True]; 
boundarydistance = 10; 
{ev, ef} = 
  NDEigensystem[{\[ScriptCapitalL], \[ScriptCapitalB]}, \[Chi][
    r], {r, 2*M, boundarydistance}, 3]; 
Print["The eigenvalues are = ", ev]
Print["The real part is = ", Re[Sqrt[ev]]]
Print["The imaginary part is = ", Im[Sqrt[ev]]]

Plot[Evaluate[ef], {r, 0, boundarydistance}]

THE PROBLEM:It works perfectly, for a condition which is not at infinity. When I try to set boundarydistance = :inf:, it gives me error.

2) NDSolve with "Shooting Method":

f = 1 - 2*(M/r); 
V = f*(l*((l - 1)/r^2) + (2*(1 - S^2)*M)/r^3); 
M = 1; l = 2; S = 2; 
\[ScriptCapitalL] = (-f)*D[f*D[\[Chi][r], r], r] + V*\[Chi][r]; 
NDSolve[{\[ScriptCapitalL] == \[Lambda]*\[Chi][r], \[Chi][2*M] == 0, \[Chi][Infinity] == 0}, \[Chi][r], {r, 0, 10}, 
  Method -> {"Shooting"}]

Plot[%, {r, 0, 10}]

THE PROBLEM: It gives the following error: NDSolve::ndsv: Cannot find starting value for the variable [Chi]^[Prime]. Besides this, how to I recover the eigenvalue? It seems like, to use NDSolve, I need to put a value for it.

3) Using the magical package in one of the answers in here:

Needs["PacletManager`"]
PacletInstall["CompoundMatrixMethod", 
   "Site" -> "http://raw.githubusercontent.com/paclets/Repository/master"]; 
Needs["CompoundMatrixMethod`"]
f = 1 - 2*(M/r); 
V = f*(l*((l - 1)/r^2) + (2*(1 - S^2)*M)/r^3); 

M = 1; l = 2; S = 2; 
sys = ToMatrixSystem[f*D[f*D[\[Chi][r], r], r] - V*\[Chi][r] == \[Omega]*\[Chi][r], 
    {\[Chi][3*M] == 0, \[Chi][Infinity] == 0}, \[Chi], {r, 3, 10}, \[Omega]]; 
Plot[Evans[\[Omega], sys], {\[Omega], 0, 15}]

THE PROBLEM: Well, it actually runs. But I need to start from somewhere far from zero (due to the singularity) and I cannot really understand the answer which is:

{{{0, 1}, {(12 - 10 r + 2 r^2 + \[FormalLambda] r^4)/((-2 + r)^2 r^2), (4 r - 2 r^2)/((-2 + r)^2 r^2)}}, {}, {}, {r, 2.5, 10}}

and a plot, which does not resembles me anything that I would like to get.

EDIT:

Since I am not converging to an answer, it may help to add some more details. The real equation that I am trying to solve is $$\frac{d^2 \psi(r)}{dr_*^2} + (\omega^2 - V(r)) \psi(r) =0$$ where the (tortoise) coordinate is defined via $$\frac{dr}{dr_*} = f(r)$$ with the definitions of both $V(r)$ and $f(r)$ are given above. If you use these two equations together, you will get in the first equation of this page.

I thought that working with a single differential equation, although a bit more complicate, might be easier. However, I have no idea, even so, how I would work with the two differential equations separately, since the problem with this tortoise coordinate is that it cannot give me a function $r(r_*)$, because it is a transcendental equation.