In case you just want to see the formula, I tried to prove, but couldn't: For $n\in \mathbb{N}$ and $n\geq 2$
$$ (n+1)^{n-1} = \frac{1}{n} \sum_{k=0}^{n-1} \binom{n}{k} (k+1)^k (n-k)^{n-1-k} $$
If you want to know, how I came to this, just read on...
My task is to derive the taylor expansion of the Lambert W function at $z=0$. My idea so far is to differentiate the equation, which defines the function $$ z=W(z)e^{W(z)}. $$ Doing so, I get that $W(z)$ satisfies the differential equation $$ W = zW'(1+W). $$ It is obvious, that $W(0)=0$, so I get $W'(0)=1.$ Now I now assume $W(z)=\sum_{k=0}^{\infty} a_n z^n$ for some radius of convergence. If $z$ is within the assumed radius, I may write this as
$$ \sum_{n=0}^{\infty} a_n z^n = z \left( \sum_{n=0}^{\infty} n a_nz^{n-1} \right) \left( \sum_{n=0}^{\infty} a_n z^n + 1\right)=\left( \sum_{n=0}^{\infty} n a_nz^{n} \right) \left( \sum_{n=0}^{\infty} b_n z^n \right) $$ with $b_k = a_k$ for $k>1$ and $b_0 = a_0 + 1 = a_0$. Using the Cauchy-Product formula I get, that the coefficients need to satisfy the condition (for $n>1$ because $a_0=0$ and $a_1 = 1$ are already known)
$$ a_n = \sum_{k=0}^{n} k a_k b_{n-k} \quad \Rightarrow \quad a_n = n a_n + \sum_{k=1}^{n-1} k a_k a_{n-k} \quad \Rightarrow \quad a_n = \frac{1}{1-n} \sum_{k=1}^{n-1} k a_k a_{n-k}$$ Applying this formula gives me to correct first terms
$$ W(z)= z -z^2 + \frac{3}{2}z^3 - \frac{8}{3}z^4 + \frac{125}{24} - \dots $$
My goal would now be, to show that indeed $a_n = \frac{(-n)^{n-1}}{n!}$. I tried to proof this by induction. So I assumed that $a_n$ is of the given form for integrers $\leq n$. I now put this into the equation from above for $a_{n+1}$ and with a little bit of elementar manipulation and a index shift it follows, that this is true, if and only if the first formula holds.