Find N degree polynomial from N+1 points

I remember in school learning how to tell if a set of numbers belonged to a quadratic function or not.

Number | first Difference | second difference
1  
         2  
3                           2  
         4  
7                           2  
         6  
13                          2  
         8  
21  

Since the second difference column was all the same, the set of numbers was said to be of 2nd order. This same method worked with third, and fourth etc. differences as well.

This led me to the conclusion that for any set of numbers, you will eventually get to a difference column of a single number. Which by definition means they're all the same.

Does this mean that given any set of n+1 points, there exists a polynomial of degree n that passes through all of them? Or does this only work if the points have a uniformly increasing x value (in this case 1). If so, how can I construct such a polynomial?

It does! And you can even make the stronger statement that the resulting polynomial is degree at most $n$ (just like with quadratics).

If you've seen systems of linear equations, there's an easy way to see why this is the case. Suppose you want a polynomial $f(x) = a_0 + a_1 x + \ldots + a_n x^n$ such that $f(x_0) = y_0$, $f(x_1) = y_1$, \ldots, $f(x_n) = y_n$ where we assume $x_i \neq x_j$ if $i \neq j$. The equations $f(x_i) = y_i$ are really $n+1$ equations in $n+1$ unknown variables $a_0,\ldots,a_n$ (i.e., in the coefficients of the unknown polynomial $f(x)$):

\begin{align} a_0 + x_0 a_1 + \ldots + x_0^n a_n &= y_0 \\ &\vdots\\ a_0 + x_n a_1 + \ldots + x_n^n a_n &= y_n \end{align}

or in matrix form

$$ \begin{pmatrix} 1 & x_0 & \cdots & x_0^n\\ 1 & x_1 & \cdots & x_1^n\\ \vdots&&&\vdots\\ 1 & x_n & \cdots & x_n^n \end{pmatrix} \begin{pmatrix} a_0\\a_1\\\vdots\\a_n \end{pmatrix} = \begin{pmatrix} y_0\\y_1\\\vdots\\y_n\end{pmatrix}. $$

The $(n+1)\times(n+1)$ matrix on the left is called the Vandermonde matrix of $x_0,\ldots,x_n$, and (with some work) you can show that its determinant is always equal to

$$ \prod_{0 \leq i < j \leq n}(x_j - x_i), $$

which is nonzero if $x_j \neq x_i$ for $j \neq i$ like we assumed for obvious reasons above. It follows that the matrix is invertible, and so there is a solution to this system of equations:

$$ \begin{pmatrix} a_0\\a_1\\\vdots\\a_n \end{pmatrix} = \begin{pmatrix} 1 & x_0 & \cdots & x_0^n\\ 1 & x_1 & \cdots & x_1^n\\ \vdots&&&\vdots\\ 1 & x_n & \cdots & x_n^n \end{pmatrix}^{-1}\begin{pmatrix} y_0\\y_1\\\vdots\\y_n\end{pmatrix}. $$

What's more, notice that we've actually proven something stronger than what we set out to. Not only is there such a polynomial of degree $n$ or less, the polynomial is unique! That is,

For any set of $n+1$ points $(x_i,y_i)$ with $x_i \neq x_j$ for $i \neq j$, there is a unique polynomial $f$ of degree at most $n$ such that $f(x_i) = y_i$ for all $i$.

If you're interested, the Wikipedia page on polynomial interpolation algorithms provides some greater technical detail about the what these coefficients look like when written out explicitly and also how efficiently some related algorithms can be implemented. This includes so-called Lagrange interpolation method, which is derived from the above solution.