I have a charge density $\rho(\vec{r})$ which is given as an expansion of spherical harmonics $Y_{l}^{m}(\hat{r})$.
$$\rho(\vec{r}) = \sum_l^{l_{max}} \sum_{m = -l}^{l} \rho_{l}(r) Y_{l}^{m}(\hat{r}) $$
Now I want to solve the Poisson equation $\nabla^2V(\vec{r}) = -4 \pi \rho(\vec{r})$. One way of solving this equation would be to use a Green's function. However, I want to explicitly solve the differential equation for each ${l,m}$.
For that, I am expanding the Poisson operator into spherical harmonics.
$$\nabla^2V(\vec{r}) = \sum_l^{l_{max}} \sum_{m = -l}^{l} \nabla^2 \left( V_l(r) \, Y_{l}^{m}(\hat{r}) \right) = \sum_l^{l_{max}} \sum_{m = -l}^{l} V_l(r) \, \nabla^2 Y_{l}^{m}(\hat{r}) + \nabla^2 V_l(r) \, Y_{l}^{m}(\hat{r}) $$
$$ \nabla^2 Y_{l}^{m}(\hat{r}) = - \frac{l (l + 1)}{r^2} \, Y_{l}^{m}(\hat{r}) $$
After some algebra, I am left with:
$$\nabla^2V_l(r) = \frac{d^2V}{dr^2} + \frac{2}{r} \frac{dV}{dr} = V_l(r) \, \frac{l (l + 1)}{r^2} - 4 \pi \rho(r)$$
This Euler-Cauchy ODE equation can than be solved. My question is now how to I get the correct initial values for solving this ODEs?