$f(z)= \frac{z}{1-e^{-z}}$ $\Leftrightarrow$ $\forall n \ge 0$ and $n\in\mathbb{Z}$, coefficient of $z^n$ in $f^{n+1}(z)$ is $1$. How to prove that if $f(z)$ is analytic in the region around origin, then $f(z)= \frac{z}{1-e^{-z}}$ $\Leftrightarrow$ $\forall n \ge 0$ and $n\in\mathbb{Z}$, coefficient of $z^n$ in $f^{n+1}(z)$ is $1$. 
I try to use the expansion of $f(z)= \sum_{n=0}^{\infty}a_n z^n$, and calculate the coefficient of $z^n$ in $(\sum_{m=0}^{\infty}a_n z^m)^{n+1}$. By this way, we can explicitly compute the first few $a_n$.
$n=0 \Rightarrow a_0=1$ 
$n=1 \Rightarrow 2 a_1 a_0 =1 \Rightarrow a_1=1/2$
$n=2 \Rightarrow 3 a_2 a_0^2 + 3 a_1^2 a_0  =1 \Rightarrow a_2=1/12$
$\cdots$
However it seems to involve a complicated partition and combination problem when $n$ is large.
That is,  firstly try to find the partition,  $\forall i, x_i\in \mathbb{Z}$ and $x_i\ge0$
$$x_1+x_2+\cdots +x_{n+1}=n\tag{1}$$
then solve the iteration function
$$\sum_{\{x_i\}} \prod_{i=1}^{n+1} a_{x_i}=1\tag{2}$$ 
with $\sum_{\{x_i\}}$ means sum over all configuration $\{x_1,\cdots, x_{n+1}\}$ such that $(1)$ holds.
How to prove the solution of $a_n$ relates to Bernoulli numbers(i.e $a_n=(-1)^n B_n/n!$ since in this question, it proves that $f(-z)=\frac{z}{e^z-1}=\sum_{n=0}^\infty B_n z^n/{n!}$)? Or by other method to prove above claim?
 A: HINT:
I think it has to do with Lagrange inversion formula for power series.
If $f= \sum_{n\ge 0} f_n z^n\in K[[z]]$ is a formal power series with $f_0=0$, $f_1 \ne 0$, then $f$ has a composition inverse $g$ such that $g(f(z)) = z$. The coefficients of $g$ can be calculated by the formula
$$g_k =\frac{1}{k} \operatorname{res} \left(\frac{1}{f}\right)^k$$
for all $k\ge 1$, where $ \operatorname{res}$ is the coefficient of $\frac{1}{z}$.
The equivalent equation you wrote is : the residue of $(\frac{1}{1-e^{-z}})^{k}$ is $1$ for all $k\ge 1$. Now you take $f(z)= 1-e^{-z}$ and calculate its formal inverse.  The inverse is $g(z) = \log \frac{1}{1-z}$, whose power series is
$$g(z) = \sum_{k \ge 1}\frac{z^k}{k}$$
A: The coefficient $[z^n]$ can be computed by integrating $f(z)/z^{n+1}$ along a closed contour winding about $z=0$. For simplicity let the contour be a circle $|z|=\epsilon$ with $\epsilon$ being small enough so that only the pole $z=0$ is inside the contour: 
$$
[z^n]=\frac{1}{2\pi i}\oint\limits_{|z|=\epsilon}\frac{dz}{(1-e^{-z})^{n+1}}
\stackrel{u=1-e^{-z}}{=}\frac{1}{2\pi i}\oint\limits_{C_\epsilon}
\frac{du}{(1-u)u^{n+1}}=1,
$$
as the residue of the integrated function in the last integral at $u=0$ is 1, which is obvious from the representation
$$
\frac{1}{(1-u)u^{n+1}}=\frac{1}{u^{n+1}}\sum_{i=0}^\infty u^i,
$$
valid for $|u|<1$. Thus the claim is proved.
The contour $C_\epsilon$ is the image of $z=\epsilon e^{i\phi}$:
$$u(\phi)=1-e^{-\epsilon e^{i\phi}}=1-e^{-\epsilon\cos\phi}[\cos(\epsilon\sin\phi)-i\sin(\epsilon\sin\phi)],
$$
which is, provided that $\epsilon<\pi$, a simple closed curve winding counterclockwise about $u=0$, with the distance to the origin varying from $1-e^{-\epsilon}$ at $\phi=0$ to $e^{\epsilon}-1$ at $\phi=\pi$. Thus, the point $u=1$ is outside the contour.
