Method 2
$\def\nn{\mathbb{N}}$
$\def\zz{\mathbb{Z}}$
$n^2$ is a degree-$2$ polynomial on $n$, so we write down the sequence $(n^2)_{n\in\nn}$ and find the differences between consecutive terms, which give rise to another sequence. Repeating this $3$ times will produce a sequence of only zeros because $3 > 2$, as we shall prove later.
0,1,4,9,...
1,3,5,...
2,2,...
0,...
Then we look at the $0$-th term in each sequence and just multiply each with the corresponding binomial coefficient to get:
$n^2 = 0 \binom{n}{0} + 1 \binom{n}{1} + 2 \binom{n}{2}$.
Again, why this is true will be proven later.
Now if we want the sum of the first $(n-1)$ squares, we just shift the coefficients to the right as follows:
$\sum_{k=0}^{n-1} = 0 \binom{n}{1} + 1 \binom{n}{2} + 2 \binom{n}{3} = \frac{1}{6}(n-1)n(2n-1)$.
Such a wonderful property of binomial coefficients is not accidental, as we shall see below.
Proof
Define the forward difference operator:
$D = ( \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$, $D(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 the following shows what happens when 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,...
Then we have:
$D\left( \text{int $n$} \mapsto \binom{n}{k+1} \right) = \left( \text{int $n$} \mapsto \binom{n}{k} \right)$ for any $k \in \zz$.
This is to be expected because it follows directly from Pascal's triangle, especially if we define $\binom{n}{k}$ using the triangle.
This means that if we have any function $f$ on $\zz$ such that $f(n) = \sum_{k=0}^\infty a_k \binom{n}{k}$ for any $n \in \zz$, then we get:
$D(f)(n) = \sum_{k=0}^\infty a_{k+1} \binom{n}{k}$ for any $n \in \zz$.
From a high-level perspective, this is the discrete version of the Taylor series, and indeed for such a function we easily see that $f(n) = \sum_{k=0}^\infty D^k(f)(0) \binom{n}{k}$ for any $n \in \zz$, because $\binom{0}{0} = 1$ while $\binom{0}{k} = 0$ for any $k \in \nn^+$.
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.
Notes
This is far more efficient that my other method in a certain sense because the series using binomial coefficients is easy to manipulate and easy to compute. Strangely, it seems that people do not appreciate methods that are better...
Also, for a wide class of non-polynomial functions, we can still compute the indefinite sum without using the series, by using the discrete analogue to integration by parts, here called summation by parts.
To derive it, simply check that $D(f \times g)(n) = f(n+1) g(n+1) - f(n) g(n) = f(n+1) D(g)(n) - D(f)(n) g(n)$ and so we get the product rule:
$D(f \times g) = R(f) \times D(g) + D(f) \times g$
where $R$ is the right-shift operator defined as:
$R = ( \text{function $f$ on $\zz$} \mapsto ( \text{int $n$} \mapsto f(n+1) ) )$
Namely for any function $f$ on $\zz$ and $n \in Z$, $R(f)(n) = f(n+1)$.
For convenience we also define the summation operator:
$S = ( \text{function $f$ on $\zz$} \mapsto ( \text{int $n$} \mapsto \sum_{k=0}^{n-1} f(k) ) )$
Then we have the important property that $DS(f) = f$ for any function $f$ on $\zz$, analogous to the fundamental theorem of calculus.
Now by substituting $f$ with $S(f)$ into the product rule and taking summation on both sides we get summation by parts:
$S( f \times g ) = S(f) \times g - S( R(S(f)) \times D(g) ) + c$ for some constant function $c$ on $\zz$.
Using this we can easily compute things like $\sum_{k=1}^n 3^k k^3$ by applying it three times, each time reducing the degree of the polynomial part.