Is there any way of characterizing the idempotent elements of the ring of polynomials in $m$ variables $x,~y,~z, \ldots$ over a finite field $F_p$, modulo $x^p-x$, so that I am really considering functions of $m$ finite variables, e.g.. $F_p[x,y,z]/\langle x^p-x, y^p-y,z^p-z \rangle$.
For example, in $F_3[x]/(x^3-x)$, $x^2$ is idempotent, because $x^4=x^3*x=x^2$.
Is there any easy way of enumerating all such idempotent elements of the ring in computer algebra systems like, magma?
The ring $\mathbb{F}_p[x_1, \dots x_m]/(x_i^p - x_i)$ is the ring of polynomial functions $\mathbb{F}_p^m \to \mathbb{F}_p$. Every function $\mathbb{F}_p^m \to \mathbb{F}_p$ is representable as a unique element of this ring, meaning this ring is actually isomorphic to the full ring of functions $\mathbb{F}_p^m \to \mathbb{F}_p$. This means it has exactly $2^{p^m}$ idempotent elements, given by those functions which take the values either $0$ or $1$.
(There are a couple ways to prove this. There's an obvious homomorphism from $\mathbb{F}_p[x_1, \dots x_m] \to \mathbb{F}_p^{\mathbb{F}_p^m}$, given by evaluations. Its kernel obviously contains $x_i^p - x_i$ for each $i$. So that gives a homomorphism between these two rings. Moreover they have the same cardinality, and the construction below shows that this homomorphism is surjective. So it is an isomorphism. We can also give a more abstract argument using the Chinese remainder theorem and tensor products.)
These functions can be more explicitly described in terms of polynomials as follows. First, it's not hard to see that the unique element of this ring which takes the value $1$ at the origin and $0$ everywhere else is
$$e_0(x_1, \dots x_m) = (1 - x_1^{p-1}) \dots (1 - x_m^{p-1}).$$
This is an idempotent. By translating this idempotent around we can produce elements which take the value $1$ at any particular point $v \in \mathbb{F}_p^m$ and $0$ everywhere else, namely (writing $v = (v_1, \dots v_m)$ and $x = (x_1, \dots x_m)$)
$$e_v(x) = e_0(x - v) = (1 - (x_1 - v_1)^{p-1}) \dots (1 - (x_m - v_m)^{p-1})).$$
These are also idempotents. (Every function $\mathbb{F}_p^m \to \mathbb{F}_p$ is a linear combination of these, so this proves surjectivity as desired above.) Altogether the $e_v$ are exactly the primitive idempotents, and the set of idempotents is exactly the set of sums of the form
$$e_S = \sum_{v \in S} e_v$$
for some subset $S \subseteq \mathbb{F}_p^m$. $e_S$ is the indicator function of the subset $S$; it takes the value $1$ on $S$ and $0$ elsewhere.
Let's see how this recovers the specific idempotent you wrote down for $p = 3, m = 1$. Here $\mathbb{F}_3[x]/(x^3 - x)$ is isomorphic to the ring of functions $\mathbb{F}_3 \to \mathbb{F}_3$. There are three primitive idempotents, namely
$$\begin{eqnarray*} e_0 &=& 1 - x^2 & & &=& (1 - x)(1 + x) \\ e_1 &=& 1 - (x - 1)^2 &=& 2x - x^2 &=& x(2 - x) \\ e_2 &=& 1 - (x - 2)^2 &=& x - x^2 &=& x(1 - x). \end{eqnarray*}$$
This gives $2^3 = 8$ idempotents coming from taking sums of subsets of these. $x^2$ itself takes the values $0, 1, 1$ at $0, 1, 2$ which means it is $e_1 + e_2$, which is not hard to verify.
An alternative way to write the primitive idempotents in one variable is as $e_a = \frac{x - x^p}{x - a}$, thinking of the factorization $x^p - x = \prod_{a \in \mathbb{F}_p} (x - a)$. We need to use $x - x^p$ instead of $x^p - x$ to get the value $1$ when we substitute $x = a$, which you can check using a formal version of l'Hopital's rule. Then you can take products of these (substituting in the different variables) to get the primitive idempotents in any number of variables, so we can rewrite $e_v$ above as
$$e_v = \prod_{i=1}^m \frac{x_i - x_i^p}{x_i - v_i}.$$
This is essentially multivariate Lagrange interpolation.