All right, this answer is a shot in the dark, but here goes.
First, transform the second-order ODE into a system of two ODEs. Let
\begin{align}
\varphi_{1} &= \psi, \\
\varphi_{2} &= \dot{\psi},
\end{align}
where the dots atop functions correspond to differentiation with respect to the independent variable (in this case, $\xi$).
Then the second-order implicit ODE
\begin{align}
\ddot{\psi}(\xi) + 2\xi^{-1}\dot{\psi}(\xi) &= e^{-\psi(\xi)} \\
\psi(0) &= 0 \\
\dot{\psi}(0) &= 0
\end{align}
can be expressed as the first-order explicit ODE
\begin{align}
\dot{\varphi}_{1}(\xi) &= \varphi_{2}(\xi) \\
\dot{\varphi}_{2}(\xi) &= -2\xi^{-1}\varphi_{2}(\xi) + e^{-\varphi_{1}(\xi)} \\
\varphi_{1}(0) &= 0 \\
\varphi_{2}(0) &= 0.
\end{align}
At first, it would appear that we cannot evaluate the right-hand side of this explicit ODE system at $\xi = 0$, like a numerical integrator requires. If a solution to this system exists, then it must be differentiable. On the assumption that a solution exists, take the limit of the right-hand side as $\xi \rightarrow 0$.
First, we know that
\begin{align}
\lim_{\xi \rightarrow 0} \varphi_{2}(\xi) = 0,
\end{align}
because we've assumed that a solution exists, so $\varphi_{2}$ is differentiable, which means it must be continuous. The limit of a continuous function at a point is its value at that point, and we know the value of $\varphi_2(0)$ because it is an initial condition.
We also know that
\begin{align}
\lim_{\xi \rightarrow 0} e^{-\varphi_{1}(\xi)} = 1
\end{align}
for similar reasons; we've assumed that $\varphi_{1}$ is differentiable, so it is continuous, and $\varphi_{1}(0) = 0$ because it is an initial condition.
Finally,
\begin{align}
\lim_{\xi \rightarrow 0} \frac{-2\varphi_{2}(\xi)}{\xi} =
\lim_{\xi \rightarrow 0} -2\dot{\varphi}_{2}(\xi),
\end{align}
by using L'Hôpital's rule on the indeterminate form $0/0$.
To proceed further, we have to make another assumption: $\dot{\varphi}_{2}$ is continuous at $\xi = 0$. Then it follows that
\begin{align}
\lim_{\xi \rightarrow 0} -2\dot{\varphi}_{2}(\xi) = -2\dot{\varphi_{2}}(0).
\end{align}
Revisiting the first-order ODE, and evaluating the right-hand side at $\xi = 0$, we can see that we have:
\begin{align}
\dot{\varphi}_{1}(0) = 0 \\
\dot{\varphi}_{2}(0) = -2\dot{\varphi}_{2}(0) + 1,
\end{align}
from which it follows that $\dot{\varphi}_{2}(0) = 1/3$.
Using this analysis, you could plug in an if statement that returns these values of the right-hand side function at $\xi = 0$, which should get you past the singularity. That said, this analysis requires a couple assumptions about continuity that may or may not hold, so take the resulting solution with a grain of salt.
