Every affine variety in $\mathbb A^n$ consisting of finitely many points can be written as the zero locus of $n$ polynomials

Assume that the points are $a_k=(a_k^{1},a_k^2,...,a_k^n)$, for $k=1,2,...,M$.

We can use the following system

$$\begin{cases}0&=\prod_{k=1}^{M}(z_1-a_k^1)\\ 0&=\prod_{k=1}^{M}(z_2-a_k^2)+\\&+\sum_{j=1}^{M}\left[\frac{(-1)^j\prod_{k=1,k\neq j}^{M}(z_1-a_k^1)}{\prod_{k=1,k\neq j}^{M}(a_j-a_k)}\cdot(z_2-a_j^2)\cdot\left\{\prod_{k=1,k\neq j}^{M}(z_2-a_k^2)+1\right\}\right]\\ ...\\ 0&=\prod_{k=1}^{M}(z_n-a_k^n)+\\&+\sum_{j=1}^{M}\left[\frac{(-1)^j\prod_{k=1,k\neq j}^{M}(z_1-a_k^1)}{\prod_{k=1,k\neq j}^{M}(a_j-a_k)}\cdot(z_n-a_j^n)\cdot\left\{\prod_{k=1,k\neq j}^{M}(z_n-a_k^n)+1\right\}\right]\end{cases}$$

The first polynomial forces the possible values for $z_1$ as $a_1^1,a_2^1,...,a_N^1$. The role of the other polynomials is to force the values of the other variables according to the value of $z_1$.

The equations are symmetric by permutations on the index $k$. Assume without loss of generality that $z_1$ is, say $=a_1^1$. Then the $r$-th equation, for $r=2,3,...,n$, becomes

$$\begin{align}0&=\prod_{k=1}^{M}(z_r-a_k^r)-(z_r-a_1^r)\cdot\left\{\prod_{k=1,k\neq j}^{M}(z_r-a_k^r)+1\right\}\\&=(z_r-a_1^r)\end{align}$$

from where $z_r$ is forced to be $=a_1^r$.

We use induction on $n$, the base case $n=1$ being trivial (note that the result is actually false for $n=0$; alternatively you can require the varieties to be nonempty and take $n=0$ as the base case).

Now suppose the result is known for $n$ and let $V\subseteq \mathbb{A}^{n+1}$ be finite. Let $a_1,\dots,a_m$ be all the different first coordinates of points of $V$ and let $V_i=\{b\in \mathbb{A}^n:(a_i,b)\in V\}$. By the induction hypothesis, for each $i$ we can choose $n$ polynomials $f_{i1},\dots,f_{in}$ whose vanishing set is $V_i$. For $1\leq k\leq n$, we can then choose a polynomial $g_k$ in $n+1$ variables such that $g_k(a_i,y)=f_{ik}(y)$ for each $i$ (here $y$ is an $n$-tuple of variables). Explicitly, if $e_i(x)$ is a polynomial in one variable that is $1$ on $a_i$ and $0$ on $a_j$ for $j\neq i$, then you can take $g_k(x,y)=\sum_i e_i(x)f_{ik}(y)$. Finally, we see that $V$ is the vanishing set of the polynomials $g_1(x,y),\dots,g_n(x,y)$ together with one more $(x-a_1)\dots(x-a_m)$.