I'm trying to found a formula for the coefficients of the power series $$ \sum_{k=0}^\infty p_k z^k = \left(\sum_{k=0}^\infty a_k z^k\right)^n. $$
Wolfram provides a nice recurrence equation for computing the set $\{p_k\}$, please see image below.
There is no reference for this, so I'm a bit skeptical about whether this is true, especially since I could not find this formula anywhere else.
The formula is correct. To prove it, let $A(z)=\sum_{n\ge 0} a_n z^n$, and let $P(z)=A(z)^n=\sum_{n\ge0}p_nx^n$. Using the chain rule, $$ P(z)' = n\big(A(z)\big)^{n-1}A'(z), $$ which after multiplying by $zA(z)$ further implies $$ A(z)\cdot zP'(z) = n\cdot zA'(z)\cdot P(z) $$ Now, using the fact that $zP'(z)=\sum_{n\ge 0}np_n z^n$ and $zA'(z)=\sum_{n\ge 0}na_n z^n$, along with the convolution formula for the product of two power series, we get that the coefficient of $z^k$ of both sides of that last equation is $$ \sum_{j=0}^k a_j \cdot (k-j)p_{k-j}=n\sum_{j=0}^kj a_{j} p_{k-j} $$ Isolating the $j=0$ term on the left, we get $$ a_0kp_k + \sum_{j=1}^k a_j \cdot (k-j)p_{k-j}=n\sum_{j=1}^kj a_{j} p_{k-j}, $$ which allows you to deduce the claimed formula.
As a side note, the recursive formula given allows you to compute $p_k$ in $O(k^2)$ time. It is also possible to do this in $O(k \log k \log n)$ time. Namely, use exponentiation by squaring to compute $\Big(\sum_{i=0}^k a_i x^i\Big)^n$ using $O(\log n)$ polynomial multiplications, only keeping track of the first $k+1$ coefficients. Each multiplication can be done in $O(k\log k)$ time using the fast Fourier transform. Wolfram's method is better when $n$ is large.