A formal power series such that $f(f(x))=x$ Let
$$f(x)=\sum_{n=1}^\infty a_n x^n$$
be a formal power series with no constant term such that $f(f(x))=x$. We find that
$$f(f(x))=a_1^2x+(a_1a_2+a_1^2a_2)x^2+(a_1a_3+2a_1a_2^2+a_3a_3)x^3+\dots$$
so $a_1^2=1$. If $a_1=1$, you need all the other terms to be zero, however if $a_1=-1$ we get a family of nontrivial solutions. Let $a_2=a$, and requiring the higher coefficients of $f(f(x))$ to be zero we can find $a_3=2a^2$, $a_4=-\frac{11}2a^3$, $a_5=\frac{11}2a^4$, $a_6=\frac{105}4a^5$, $a_7=-\frac{279}2 a^6$...
Is there a closed form for these numbers?

The claim that the coefficients are unique given the quadratic term is incorrect. Before going into that, here's some background on Anne's answer. We consider the simpler question of how to determine solutions to $f(f(x)) = x$ where $f$ is a Mobius transformation $f(x) = \frac{ax + b}{cx + d}$. Over a field $K$ the group of Mobius transformations is isomorphic to the projective general linear group $PGL_2(K)$, with the isomorphism given by sending a $2 \times 2$ matrix $\left[ \begin{array}{cc} a & b \\ c & d \end{array} \right]$ to $\frac{ax + b}{cx + d}$. So the problem reduces to finding matrices $M$ squaring to a scalar multiple of the identity.
So write $M^2 = c$ for some scalar $c$, so that the eigenvalues of $M$ are $\pm \sqrt{c}$ (over an algebraically closed field, say $\mathbb{C}$). Since we only need to work up to scale we might as well divide $M$ by $c$, or equivalently assume WLOG that $c = 1$, so $M^2 = 1$ and $M$ has eigenvalues $\pm 1$. If $M$ is conjugate to a nontrivial Jordan block then it squares to another nontrivial Jordan block so can't square to $1$; hence $M$ is diagonalizable. If $M$ has eigenvalues $\{ 1, 1 \}$ or $\{ -1, -1 \}$ then it's a scalar multiple of the identity; otherwise its eigenvalues are $\{ 1, -1 \}$. This implies that $\text{tr}(M) = 0$ and $\det(M) = -1$, so $M$ has the form $\left[ \begin{array}{cc} a & b \\ c & -a \end{array} \right]$ where $\det(M) = -a^2 - bc = -1$. This gives $bc = 1 - a^2$.
Now we add the constraint that as a formal power series $f(x)$ should have no constant term. This means $f(0) = 0$ which gives $b = 0$. Then $a = \pm 1$ and we can take $a = 1$ WLOG, which gives $M = \left[ \begin{array}{cc} 1 & 0 \\ c & -1 \end{array} \right]$, hence
$$f(x) = \frac{x}{cx - 1} = -x - cx^2 - c^2 x^3 - \dots $$
as in Anne's answer (with $a = -c$).

On the other hand Klaus's answer implies that the resulting power series cannot be unique, since we can conjugate by an invertible (with respect to composition) formal power series. With a little more effort (showing that we can even arrange for the first $n$ coefficients of this series past $x$ to vanish) we can show that the first $n$ coefficients never uniquely determine the others. So you did something funny with your calculations but I'm not sure what.
Generally, it's known that every formal power series of finite order (with respect to composition) is conjugate to $f(x) = \zeta x$ for $\zeta$ some root of unity (above we have $\zeta = -1$, and in general we probably need to work over an algebraically closed field). This implies:

Classification: If $f(x)$ is a formal power series satisfying $f(0) = 0$ and $f(f(x)) = x$, then either $f(x) = x$ or $\boxed{ f(x) = g(-g^{-1}(x)) }$ where $g(x) = x + \dots$.

Given $g(x)$, the coefficients of $g^{-1}(x)$ can be computed using Lagrange inversion. This gives solutions depending on an infinite number of parameters, namely the higher coefficients $g_i$ of $g(x)$, and the first $n$ coefficients of $f$ depend only on the first $n$ coefficients of $g$. For example, if $g(x) = x + g_2 x^2 + \dots$ then $g^{-1}(x) = x - g_2 x^2 + \dots$ which gives $f(x) = -x + 2 g_2 x^2 + \dots$.
To give a relatively explicit example, take $g(x) = x - ax^2$. Then $g^{-1}(x) = \frac{1 - \sqrt{1 - 4ax}}{2a}$ by the quadratic formula (we have to take the minus sign so that $g^{-1}(0) = 0$); this is a version of the generating function of the Catalan numbers. If $y - ay^2 = x$ then $g(-y) = -y - ay^2 = x - 2y$, which gives
$$\boxed{ \begin{align*} f(x) &= x - \frac{1 - \sqrt{1 - 4ax}}{a}  \\
 &= -x - 2ax^2 - 4a^2 x^3 - 10a^3 x^4 - \dots \\
 &= -x - \sum_{n=2}^{\infty} \frac{2}{n} {2n-2 \choose n-1} a^{n-1} x^n. \end{align*} }$$
