Background
Long ago I bumped into an exercise in ordinary differential equations, which asks to find a solution to the differential equation:
$$h'(x)=\frac{1}{2(1+xh(x))}$$
It turns out that $h(x)$ is the inverse function of
$$G(x)=2e^{x^2}\int_0^xe^{-u^2}\,du=\sqrt{\pi}e^{x^2}\hbox{erf}(x)$$
where $\hbox{erf}(x)$ is the standard error function. This function $G(x)$ (multiplied by $1/2$), appears in the problem set of the second chapter of the book "Special functions and their applications", by Lebedev, where the reader is asked to prove that it satisfies a certain differential equation, namely,
$$G'(x)=2(1+xG(x))$$
whence it easily follows that the inverse $G^{-1}(x)=h(x)$ satisfies the ODE above for $h$. In Lebedev's book the reader is asked to use the differential equation for $G$ to derive a series expansion for the error function, valid in the entire complex plain. In terms of $G$, the series is as follows:
$$G(z)=\sum_{k=0}^{\infty}\frac{2^kz^{2k+1}}{(2k+1)!!}\qquad (z\in\mathbb{C})$$
Lebedev leaves it at that, but some further comments may be of interest.
A simple calculation involving Stirling's formula shows that this series has an infinite radius of convergence, hence it defines an entire function in the complex plain. Moreover, restricting to positive real values only, the series expansion shows that all derivatives of $G$ are positive in $(0,\infty)$, i.e., $G(x)$ is an absolutely monotone function in $(0,\infty)$. Moreover, since $G'(0)\neq 0$, there is a disc around zero where the inverse function $G^{-1}(z)=h(z)$ is analytic, and its series expansion can be computed from that of $G(z)$. However, it turns out that the radius of convergence of the series for $G^{-1}(z)$ around zero is finite. This perhaps is connected to the fact that $G(z)$ is not one-to-one in the whole complex plain, i.e., does not have a global inverse.
The problem
Looking into various properties of the inverse function $h(x)=G^{-1}(x)$, I came across the function $g(x)=\exp(G^{-1}(x)^2)=\exp(h^2(x))$, where $x\geq 0$.
Recall that a function $f(x)$ defined on $(0,\infty)$ is called completely monotone if it has derivatives of all order and $(-1)^nf^{(n)}(x)\geq 0$ for all $n=0,1,2,\dots$.
The problem is this:
Is the function $g''(\sqrt{x})$ completely monotone in $[0,\infty)$?
Using the differential equation for $h$ above, the second derivative is
$$g''(x)=\frac{\exp(h^2(x))}{2(1+xh(x))^3}$$ so $$g''(\sqrt{x})=\frac{\exp(h^2(\sqrt{x}))}{2(1+\sqrt{x}h(\sqrt{x}))^3}$$
Since this is a composition of $g''(x)$ and the square root function $\sqrt{x}$, and since the derivative of $\sqrt{x}$ is completely monotone, then if we knew that $g''(x)$ is itself completely monotone, we would be able to deduce that $g''(\sqrt{x})$ is completely monotone as well. But that is not the case: $g''(x)$ is not completely monotone, because, for example, it is not convex. Its second derivative, the fourth derivative of $g(x)$, is not a positive function.
A partial result
With $f(x)=g''(\sqrt{x})$, the following holds. For every $n$, there is some $x_n>0$, such that for all $x>x_n$ we have $(-1)^nf^{(n)}(x)>0$.
My proof is neither elegant nor short. Moreover, I still don't know whether $(-1)^nf^{(n)}(x)>0$ holds for all $x$.
There is lots of information and literature involving completely monotone functions. For example, the property of being completely monotone is equivalent to the property of being the Laplace transform of a positive measure, and in our case, since $g''(0)=\frac{1}{2}<\infty$, we would have a bounded positive measure. There are lots of tricks and special cases that help determine whether a function is completely monotone, but as far as my research has reached, I could not settle the case for this specific function. It might be interesting to discover new methods of proving -- or disproving -- that a certain function is completely monotone.
Any insights, comments, remarks and proofs are welcome.