Consider a continuous real-valued monotone increasing function $f:\mathbb R^+\to\mathbb R^+$ satisfying $f\!\left(x+x^\phi\right)=x,$ where $\phi=\frac{1+\sqrt5}2$ is the golden ratio. Here is a plot of the function $f(z)$:
I need to find an asymptotic expansion of $f(z)$ for $z\to0^+$ in terms of powers of $z$. By manually balancing coefficients, I was able to find a few first terms: $$f(z)=z-z^\phi\color{red}+\phi\,z^{2\,\phi-1}+\mathcal O\left(z^{3\,\phi-2}\right), \quad z\to0^+.$$ How can I find more terms and, ideally, a general formula for all terms of this series?
Update: A wrong sign in formula is corrected (red).
I notice that factoring $z$ out yields $f^{-1}(z)=z\big(1-z^{1/\varphi}\color{Red}{+}\varphi z^{2/\varphi}+\mathcal{O}(z^{3/\varphi})\big)$. Thus, let's write
$$ f^{-1}(z)=z\sum_{k=0}^{\infty} a_k z^{k/\varphi} $$
for some coefficients $a_k$ with $a_0=1$. Substitute $z=x+x^{\varphi}=x(1+x^{1/\varphi})$ to get
$$ x=x(1+x^{1/\varphi})\sum_{k=0}^{\infty} a_k(1+x^{1/\varphi})^{k/\varphi}x^{k/\varphi}. $$
Divide by $x(1+x^{1/\varphi})$, then expand Newton-Binomial series:
$$ \frac{1}{1+x^{1/\varphi}} = \sum_{k=0}^{\infty} a_kx^{k/\varphi}(1+x^{1/\varphi})^{k/\varphi} $$
$$ \begin{array}{ll} \displaystyle \sum_{n=0}^\infty (-1)^nx^{n/\varphi} & \displaystyle = \sum_{k=0}^{\infty} a_k x^{k/\varphi}\sum_{j=0}^{\infty} \binom{k/\varphi}{j}x^{j/\varphi} \\ & \displaystyle =\sum_{n=0}^\infty x^{n/\varphi} \left[\sum_{k=0}^n \binom{k/\varphi}{n-k} a_k\right]. \end{array} $$
Equating coefficients yields a recurrence for $a_n$:
$$ a_n=(-1)^n-\sum_{k=1}^{n-1} \binom{k/\varphi}{n-k}a_k. $$
This also shows there's a sign error in your quadratic term.