During my physics research in cosmology I encountered a differential equation of the following type ($f,g$ are functions that only depend on the spherical coordinate $r$, $\partial_i$ means $\partial/\partial x^i$ where $x^i \in \{x,y,z\}$ etc.):
\begin{equation} \nabla^2 f - [a_1(\nabla^2 f)^2 - a_2 \sum^3_{i,j=1}(\partial_i \partial_j f)(\partial_i \partial_j f)] = -g(r) \end{equation}
where $a_1,a_2$ are arbitrary (non-zero) real numbers (constants).
Useful identities are (since $f = f(r)$):
\begin{equation} (\nabla^2 f)^2 = \Big(\frac{d^2 f}{dr^2}\Big)^2 + \frac{4}{r}\Big(\frac{d^2 f}{dr}\Big)\Big(\frac{df}{dr}\Big) + \frac{4}{r^2} \Big(\frac{df}{dr}\Big)^2 \end{equation}
and
\begin{equation} \sum^3_{i=1} \sum_{j = 1}^3 (\partial_i \partial_j f)(\partial_i \partial_j f) = \frac{2}{r^2}\Big(\frac{df}{dr}\Big)^2 + \Big(\frac{d^2 f}{dr^2}\Big)^2 \end{equation}
We would like to get an analytical solution for for $df/dr$ if it exists and otherwise a numerical one would be useful.