I will refer to Qiaochu's excellent answer here as proof that if we define

$$f(N):=\sum\limits_{n=0}^N n^2$$

then $f$ is a polynomial of degree $3$.

It is easy to calculate the first few values of this sum. Namely,

$\begin{align} f(0) &= 0 \\ f(1) &= 1 \\ f(2) &= 5 \\ f(3) &= 14 \end{align}$

I claim that these four points are sufficient to uniquely determine $f$.

To wit, we have in general that

$$f(x)=\sum\limits_{k=0}^3 c_kx^k$$

which when be combined with the four computed values above results in the following system of equations:

$$\begin{pmatrix} 1 & 0 & 0 & 0 \\ 1 & 1 & 1 & 1 \\ 1 & 2 & 4 & 8 \\ 1 & 3 & 9 & 27 \end{pmatrix} \begin{pmatrix} c_0 \\ c_1 \\ c_2 \\ c_3 \end{pmatrix} = \begin{pmatrix} 0 \\ 1 \\ 5 \\ 14 \end{pmatrix}$$

This matrix is a Vandermonde matrix which has a well-known determinant

$$\begin{align} \det(V) &= (1-0)(2-0)(3-0)(2-1)(3-1)(3-2) \\ &\neq 0 \end{align}$$

Because its determinant is nonzero, the matrix is invertible, and so we have

$$\begin{pmatrix} c_0 \\ c_1 \\ c_2 \\ c_3 \end{pmatrix} = V^{-1}\cdot\begin{pmatrix} 0 \\ 1 \\ 5 \\ 14 \end{pmatrix}$$

from which $f(x)$ can be determined directly.

However, if you're like most people, inverting a $4\times4$ matrix doesn't exactly tickle your fancy!

Luckily, now that we see that the interpolating cubic is unique, we could find it through the described matrix multiplication, but we would get to the same result if we proceeded a different route as well. This is where Lagrange polynomials come to the rescue.

Using the general formula, we have immediately that

$$\begin{align} f(x) &= 0\cdot(\dots)+1\cdot\frac{x(x-2)(x-3)}{1(1-2)(1-3)}+5\cdot\frac{x(x-1)(x-3)}{2(2-1)(2-3)}+14\cdot\frac{x(x-1)(x-2)}{3(3-1)(3-2)} \\ &= \frac{1}{2}\left(x^3-5x^2+6x\right)-\frac{5}{2}\left(x^3-4x^2+3x\right)+\frac{14}{6}\left(x^3-3x^2+2x\right) \\ &= \frac{1}{6}\left(2x^3+3x^2+x\right) \\ &= \frac{1}{6}x\left(x+1\right)\left(2x+1\right) \end{align}$$

You can generalize this approach to find expressions for $\sum n^p\quad\forall p\in\mathbb{N}$.

Or, you know, there's always Faulhaber's formula.

I will refer to Qiaochu's excellent answer here as proof that if we define

$$f(N):=\sum\limits_{n=0}^N n^2$$

then $f$ is a polynomial of degree $3$.

It is easy to calculate the first few values of this sum. Namely,

$\begin{align} f(0) &= 0 \\ f(1) &= 1 \\ f(2) &= 5 \\ f(3) &= 14 \end{align}$

I claim that these four points are sufficient to uniquely determine $f$.

Since $$\deg f = 3$$, we have in general that

$$f(x)=\sum\limits_{k=0}^3 c_kx^k$$

which when be combined with the four computed values above results in the following system of equations:

$$\begin{pmatrix} 1 & 0 & 0 & 0 \\ 1 & 1 & 1 & 1 \\ 1 & 2 & 4 & 8 \\ 1 & 3 & 9 & 27 \end{pmatrix} \begin{pmatrix} c_0 \\ c_1 \\ c_2 \\ c_3 \end{pmatrix} = \begin{pmatrix} 0 \\ 1 \\ 5 \\ 14 \end{pmatrix}$$

Note that $c_0 = 0$ trivially, and so we can help ourselves by writing the reduced matrix equation for the other three coefficients as

$$\begin{pmatrix} 1 & 1 & 1 \\ 2 & 4 & 8 \\ 3 & 9 & 27 \end{pmatrix} \begin{pmatrix} c_1 \\ c_2 \\ c_3 \end{pmatrix} = \begin{pmatrix} 1 \\ 5 \\ 14 \end{pmatrix}$$

Call the above matrix $V$. This matrix is a Vandermonde matrix which has a well-known determinant

$$\begin{align} \det(V) &= 1\cdot 2\cdot 3\cdot(2-1)(3-1)(3-2) \\ &\neq 0 \end{align}$$

Because its determinant is nonzero, the matrix is invertible, and so we have

$$\begin{pmatrix} c_1 \\ c_2 \\ c_3 \end{pmatrix} = V^{-1}\cdot\begin{pmatrix} 1 \\ 5 \\ 14 \end{pmatrix}$$

from which $f(x)$ can be determined directly.

However, if you're like most people, inverting matrices doesn't exactly tickle your fancy!

Luckily, now that we see that the interpolating cubic is unique, we could find it through the described matrix multiplication, but we would get to the same result if we proceeded a different route as well. This is where Lagrange polynomials come to the rescue.

Using the general formula, we have immediately that

$$\begin{align} f(x) &= 0\cdot(\dots)+1\cdot\frac{x(x-2)(x-3)}{1(1-2)(1-3)}+5\cdot\frac{x(x-1)(x-3)}{2(2-1)(2-3)}+14\cdot\frac{x(x-1)(x-2)}{3(3-1)(3-2)} \\ &= \frac{1}{2}\left(x^3-5x^2+6x\right)-\frac{5}{2}\left(x^3-4x^2+3x\right)+\frac{14}{6}\left(x^3-3x^2+2x\right) \\ &= \frac{1}{6}\left(2x^3+3x^2+x\right) \\ &= \frac{1}{6}x\left(x+1\right)\left(2x+1\right) \end{align}$$

You can generalize this approach to find expressions for $\sum n^p\quad\forall p\in\mathbb{N}$.

Or, you know, there's always Faulhaber's formula.