Here is the ODE I want to integrate, $$R''(y)-\frac{2}{k-y}R'(y)-\frac{l(l+1)}{(k-y)^{2}}R(y)=0$$ We see that it has a regular singular point at $y=k$ where $k<0$. Is there a way to obtain the solution of this ODE? What I am thinking is to obtain a series solution of the ODE of the form,
$$R(y)=c_{1}(y-k)^{\alpha_{1}}(a_{1}+a_{2}(y-k)+a_{3}(y-k)^{2}+...)+c_{2}(y-k)^{\alpha_{2}}(b_{1}+b_{2}(y-k)+b_{3}(y-k)^{2}+...)$$
and then find the indicial exponents $\alpha_{1}$ and $\alpha_{2}$ as well as the constants $a_{1},a_{2},...$ and $b_{1}, b_{2},...$. The crucial part is the computation of the constants $c_{1}$ and $c_{2}$ which I think can be obtained numerically by sampling points in the neighborhood of $y=k$ and then linearly fit with $R$ to obtain $c_{1}$ and $c_{2}$. But I am not sure of that. Is there a better efficient computational method of integrating the ODE?
Plugging the standard ansatz $R(y) = (y-k)^r$ into the ODE, we get the characteristic polynomial
$$ r(r-1) + 2r - l(l+1) = 0 $$ $$ \implies (r-l)(r+l+1) = 0 $$ $$ \implies r = l, -1-l $$
So the general solution is
$$ R(y) = c_1(y-k)^l + c_2(y-k)^{-1-l} $$
The constants are obtained from the initial conditions, which you do not have.