This is a simple presentation of the Newton series that has been mentioned in *user90369*'s answer.

The forward difference operator $Δ$ defined as:
$
\def\nn{\mathbb{N}}
\def\zz{\mathbb{Z}}
\def\lfrac#1#2{{\large\frac{#1}{#2}}}
\def\lbinom#1#2{{\large\binom{#1}{#2}}}
$

$Δ = ( \text{function $f$ on $\zz$} \mapsto ( \text{int $n$} \mapsto f(n+1) - f(n) ) )$

Namely, for any function $f$ on $\zz$ and $n \in \zz$, we have $Δ(f)(n) = f(n+1) - f(n)$.

If you think of the functions as sequences (infinite in both directions), then taking the forward difference means replacing each term with the value of the next term minus itself. For example if you repeatedly take the forward difference of the sequence of cubes:

```
...,-27,-8,-1, 0, 1, 8,27,...
..., 19, 7, 1, 1, 7,19,37,...
...,-12,-6, 0, 6,12,18,24,...
..., 6, 6, 6, 6, 6, 6, 6,...
..., 0, 0, 0, 0, 0, 0, 0,...
..., 0, 0, 0, 0, 0, 0, 0,...
```

This powerful abstraction makes it easy to get a lot of things. For instance, the numbers obtained here can be easily used to obtain the general formula for sum of cubes, as you desire.

**General method for indefinite summation**

The key is that:

$Δ\left( \text{int $n$} \mapsto \lbinom{n}{k+1} \right) = \left( \text{int $n$} \mapsto \lbinom{n}{k} \right)$ for any $k \in \zz$.

This is to be expected because it follows directly from extending Pascal's triangle naturally, namely if we define $\lbinom{n}{k}$ by the recurrence:

$\lbinom{n}{0} = 1$ for any $n \in \zz$.

$\lbinom{0,k+1}{0} = 0$ for any $k \in \nn$.

$\lbinom{n+1}{k+1} = \lbinom{n}{k+1} + \lbinom{n}{k}$ for any $k \in \nn$ and $n \in \zz$.

Now consider any function $f$ on $\zz$ such that $f(n) = \sum_{k=0}^d a_k \lbinom{n}{k}$ for any $n \in \zz$. Then we have for any $m \in \nn_{\le d}$:

$Δ^m(f)(n) = \sum_{k=0}^{d-m} a_{k+m} \lbinom{n}{k}$ for any $n \in \zz$.

And hence:

$Δ^m(f)(0) = a_m$.

Which immediately gives Newton's series:

$f(n) = \sum_{k=0}^d Δ^k(f)(0) \lbinom{n}{k}$ for any $n \in \zz$.

From a high-level perspective, this is the discrete version of the Taylor series.

This works for any polynomial function $f$ on $\zz$, since $D^k(f)$ is the zero function once $k$ is larger than the degree of $f$, so we can use it to immediately find the series for $(\text{int n} \mapsto n^3)$, and then just take the anti-difference by shifting the coefficients of the series the other way. The undetermined constant that appears will drop out once we perform a definite sum like if we want the sum of the first $m$ cubes.

**Sum of $p$ powers**

For example if we want $\sum_{k=1}^{n-1} k^3$ we first find the iterated forward differences of the sequence of cubes $( n^3 )_{n \in \zz}$:

```
..., 0, 1, 8,27,64,...
..., 1, 7,19,37,...
..., 6,12,18,...
..., 6, 6,...
..., 0,...
```

So we immediately get $n^3 = 0 \binom{n}{0} + 1 \binom{n}{1} + 6 \binom{n}{2} + 6 \binom{n}{3}$ and hence $\sum_{k=0}^{n-1} k^3 = 0 \binom{n}{1} + 1 \binom{n}{2} + 6 \binom{n}{3} + 6 \binom{n}{4} = \lfrac{n(n-1)}{2} \Big( 1 + \lfrac{6(n-2)}{3} \big( 1 + \lfrac{n-3}{4} \big) \Big) = \Big( \lfrac{n(n-1)}{2} \Big)^2$.

**Computation efficiency**

This is far more efficient than the typical method in some textbooks (namely by taking summation on both sides of $(n+1)^3-n^3 = 3n^2+3n+1$ and telescoping) because the Newton series is easy to compute and manipulate. For sum of $p$-powers we only need $O(p^2)$ arithmetic operations to find the forward-differences and then $O(p^2)$ more to simplify the series form into a standard polynomial form. In contrast, the other method requires $O(p^3)$ arithmetic operations.