Understanding the singularity in $f^{-1}(x)=\int_0^x f(t)dt$

Note that here $f^{-1}(x)$ is the functional inverse of $f$. Clearly from the definition of the equation $f^{-1}(0)=0\Rightarrow f(0)=0$. By setting $x\to f(x)$ and differentiating one finds

$$f'(x)f(f(x))=1$$

and therefore $\lim_{x\to 0}f'(x)=\infty$. Given this peculiar behavior I decided to try an ansatz $f(x)=Cx^a$ which does yield a solution upon substitution with $a=\phi-1$ , $C=(a+1)^{a/(a+1)}$ with $\phi$ being the golden ratio.

I was wondering if this solution is unique, or can we find another solution satisfying the same initial condition? It seems to me that $f(0)=0$ is a very peculiar initial condition that could potentially give rise to non-analytic solution(s) around the origin. In the differentiated form of the equation, the problem does not arise if $f(0)\neq 0$, in which case a regular series solution can be developed. Is there a name for this phenomenon of singular initial conditions in the context of nonlinear ODEs?

This is a step towards a series solution for a subset of solutions ($f(x)$ such that $f(x)=x$ has a solution for $x \not = 0$). Let $f(x_0)=x_0$ for some $x_0 \not = 0$. Then you can find a series solution centered around $x_0$.

From $$f'(x)f(f(x))=1 \tag 1$$

you can find that $f'(x_0)=\frac{1}{x_0}$.

Differentiating $(1)$ yields $$f'(x)f(f(x))=1$$ yields $$f''(x)f(f(x))+(f'(x))^2f'(f(x))=0$$

From which you can get that $f''(x_0)=-\frac{1}{x_0^4}$. Letting $a_n=f^{(n)}(x_0)$, a recursive formula can be obtained from the general Leibniz rule and Faa di Bruno's formula by differentiating $(1)$ a total of $n-1$ times. Altogether, this gives the "result" (result in quotes because this is so unwieldy) of $$a_n=-\frac{1}{x_0}\sum_{k=1}^{n-1}\binom{n-1}{k}a_{n-k}\sum_{i=1}^ka_i B_{k,i}(a_1,\cdots, a_{k-i+1})$$

where $B_{k,i}$ represents the Bell polynomial, and $a_1=\frac{1}{x_0}$. Some initial values are $a_2=-\frac{1}{x_0^4}$, $a_3=\frac{1+3x_0}{x_0^8}$, $a_4=-\frac{1+3x_0+10x_0^2+15x_0^3}{x_0^{13}}$, $a_5=\frac{1+3x_0+10x_0^2+30x_0^3+55x_0^4+105x_0^5+105x_0^6}{x_0^{19}}$. Then $$f(x)=x_0+\sum_{n=1}^{\infty}\frac{a_n}{n!}(x-x_0)^n$$

and plugging in $x = 0$ should yield $0$.

Graphically, it seems like the only solution is with $x_0=\phi$, which would give the solution you got of $f(x)=\phi^{(\phi-1)/\phi}x^{\phi-1}$. This would mean (if true) that all other solutions don't cross the line $y=x$, except at $x=0$.

Note that if $f(x)$ is continuous (and differentiable) everywhere, then it must be monotonic for the inverse to exist. If $f(x)$ is decreasing everywhere, then $f(x) > 0$ for $x<0$ and $f(x)<0$ for $x>0$ since $f(0) = 0$. But $f'(x)f(f(x))<0$ for $x>0$, so it can't equal $1$. If $f(x)$ is increasing everywhere, then $f(x) < 0$ for $x<0$ and $f(x)>0$ for $x>0$. But $f'(x)f(f(x)) < 0$ for $x<0$.

Thus, there's no everywhere differentiable $f(x)$ that satisfies the relation.

If $f(x)$ is just defined on $x \ge 0$ (and differentiable for $x>0$), then there's three cases: it intersects $y = x$ at $x=\phi$ (which leads to your $f(x)$), $0 <f(x)<x$ for all $x>0$, or $f(x)>x$ for all $x > 0$. If $f(x)>x$, then $f(f(x))>x$ and $0<f'(x) < \frac{1}{x}$. But this means that $\lim_{x \to \infty}f'(x) = 0$, which is impossible to have with $f(x)>x$.

If $0<f(x)<x$, $f(f(x))<x$ and $\frac{1}{x} < f'(x)$. But this means that $f'(x) \to \infty$ as $x \to 0$, which is impossible to have with $0<f(x) < x$.

Thus, your solution is the only solution that's differentiable and continuous for $x>0$.