Why does this trick gives precisely the formulas for the sum of the $n$ first naturals and the $n$ first squares?

Suppose that $f(x)$ is a polynomial of degree $n$. We denote $\Delta f(x)$ to represent $f(x+1)-f(x)$, and call this the forward difference operator. Using the binomial theorem we have $\Delta f(x)$ is of degree $n-1$. Using the binomial theorem again we deduce that $\Delta \Delta f(x)$ is of degree $n-2$.......$\Delta^{n-1}f(x)$ is of degree $1$, hence $\Delta^{k} f(x)=0$ for $k > n$.

Now suppose we have the sequence $f(0),f(1),f(2),...f(n)$ and we want to find an explicit form for $f(x)$. Then, note that we can sum both sides of:

$$f(i+1)-f(i)=\Delta f(i)$$

From $1$ to $x-1$ to get $f(x)=f(0)+\sum_{i=0}^{x-1} \Delta f(i)$ as we have a telescoping series. Now we can represent $\Delta f(x)$ in a similar manner, $\Delta \Delta $ in a similar manner ..... and we can combine all this.From this we can see that if the series converges:

$$f(x)=f(0)+\Delta (0) \sum_{x_0=0}^{x-1} 1+\Delta \Delta (0) \sum_{x_0=0}^{x-1} \sum_{x_1=0}^{x_0-1} 1+\Delta \Delta \Delta (0) \sum_{x_0=0}^{x-1} \sum_{x_1=0}^{x_0-1} \sum_{x_2=0}^{x_1-1} 1+\cdots$$

Where $\Delta^i (0)$ denotes the first term ($x=0$) of the $i$ th difference sequence of $f(x)$.

Through a combinational argument, if we take ${x \choose 0}=1$ and $\Delta^{0}(0)=f(0)$ one can show that is equal to:

$$\sum_{i=0}^{\infty} \Delta^i(0){x \choose i}$$

If we assume $\Delta^n(f(x))=c$, which is a safe assumption given the $n+1$ data points (take the difference sequences and see what I mean, it leads to a correct answer but not all correct answers ), we have that $\Delta^i(0)=0$ for $i >n$ so this equation reduces to:

$$\sum_{i=0}^{n} \Delta^i(0){x \choose i}$$

The important thing to note is that the leading coefficient is of most degree $n$.

$\Delta^i(0)=0$ for $i>n$ follows because one can assume $\Delta^n(f(x))=c$ given these $n$ terms (take the finite differences and see what I mean), but what if $c=0$? Then it can be deduced that $f(x)$ is degree at most $n-1$. On the other hand if $c \neq 0$ then from the equation we must have that $f(x)$ is of degree $n$.Furthermore if $c=0$, then we also have that $\Delta^{n-1}(f(x))=c_2$. And if $c_2=0$ then $f$ is at most degree $n-2$. So the important thing is that if $f$ is at minimum $d$ then $d$ is the first positive integer such that $\Delta^d(f(x))=c \neq 0$. Furthermore if $\Delta^n(f(n))=c \neq 0$ then we must have that $f$ is of degree $n$.

Note that:

$$\Delta^3(s(x))=\Delta (\Delta (\Delta \sum_{i=1}^{x} i^2)))=\Delta (\Delta (\Delta s(x)))=\Delta \Delta (x+1)^2=\Delta (2x+3)=2 \neq 0$$

We must have then that $s$ is of degree $3$. Adding values won't change a thing so long as they follow $s(x)$ because $\Delta^i(0)=0$ for $i>3=\text{degree}$ (note I haven't assumed that$\Delta^3(s(x))=c \neq 0$, rather I proved it). On the other hand if we add values that do not follow $s(x)$ then we may find a different (higher degree) polynomial (again think of the differences). Thus your polynomial is of minimum degree for $4$ given terms (it can be of larger degree just take the next value in such a way that it does not follow $s(x)$), and we call this a Lagrange polynomial. Let me also note that given $n+1$ data points $f(0),f(1),...f(n)$ if we assume $\Delta^n (f(x))=c$ then we consequently find the Lagrange polynomial that corresponds to this sequence.

A line is described uniquely using $2$ points.

A parabola is described uniquely using $3$ points.

A cubic is described uniquely using $4$ points.

In general, a polynomial with degree $n$ is described uniquely using exactly $n+1$ points.

Since you give the first $4$ points, it will uniquely describe the cubic polynomial that produces the sum of the first $n$ squares. You will notice that if you remove any point from the cubic or reduce the degree of the polynomial, it will no longer produce the same function.

This is an example of working out Ahmed's answer, which is good and contains a lot of information.

This trick works best if you start from $0$. To find the polynomial $f(x)$ with $$f(0) = 0, \; f(1) = 1, \; f(2) = 5, \; f(3) = 14, \; f(4) = 30,$$ form a triangle of differences

\begin{align*} &0\quad \quad 1 \quad\quad 5 \quad\quad 14 \quad\quad 30 \\ &\quad 1 \quad\quad 4 \quad\quad 9 \quad\quad 16 \\ &\quad\quad 3 \quad\quad 5 \quad\quad 7 \\ &\quad\quad\quad 2 \quad\quad 2 \\ &\quad\quad \quad\quad 0\end{align*}

and take the left side of the triangle as coefficients in a Newton series:

\begin{align*} f(x) &= 0 \binom{x}{0} + 1 \binom{x}{1} + 3 \binom{x}{2} + 2 \binom{x}{3} + 0 \binom{x}{4} \\ &= 0 + x + \frac{3x (x-1)}{2} + \frac{2 x(x-1)(x-2)}{6} \\ &= \frac{x}{6} + \frac{x^2}{2} + \frac{x^3}{3}.\end{align*}