I know that for $$\left(\sum_{n=0}^{\infty} a_n x^n\right)^2=\sum_{n=0}^{\infty} \left(\sum_{j=0}^{n} a_j a_{n-j}\right)x^n,$$ however I'm not sure precisely how to extend this to a cube rather than a square. Certainly it'll be applying the above identity again, but I think the double sum is tripping me up!
I would really appreciate some help showing me how to index my sums at the end here!
If we rewrite the formula for squaring in a slightly more symmetric way, it becomes slightly easier to see the pattern:
$$\left(\sum_{n=0}^{\infty} a_n x^n\right)^2 = \sum_{n=0}^{\infty} \left(\sum_{\substack{i+j=n \\ i,j \geq 0}} a_i a_j \right)x^n.$$
You can work out a similar formula for cubing, or taking 4th powers, etc.
$$\left(\sum_{n=0}^{\infty} a_n x^n\right)^3 = \sum_{n=0}^{\infty} \left(\sum_{\substack{i+j+k=n \\ i,j,k \geq 0}} a_i a_j a_k \right)x^n.$$
If you don't like the idea of summing over a subset of ordered triples, one can rewrite things in the following way, which is perhaps easier for a computer to implement but which I personally find more difficult to understand at a glance.
$$ \sum_{\substack{i+j+k=n \\ i,j,k \geq 0}} a_i a_j a_k = \sum_{i=0}^n \sum_{j=0}^{n-i} a_i a_j a_{n-i-j}. $$
The idea behind the conversion is that a priori we only know that i is between 0 and n, but for a given value of i, that puts more restrictions on j, and once we know both i and j, we know what k will be. This is completely analogous to how one expresses an integral over a region as a double integral or triple integral. With both sums and integrals, it is often conceptually convenient to express things over the region or set until actual computations need to be done, and only then to expand into multiple sums or multiple integrals.