Does this hold for $n\times n$ matrices?
Of course not. The problem statement seems to insinuate that $A$ and $B$ may possibly not commute when one of them has a zero trace. So, suppose you can find an example in which both $A$ and $B$ have zero trace and $AB\ne BA$. Then $\pmatrix{1\\ &A}$ and $\pmatrix{1\\ &B}$ will serve as a counterexample for $3\times3$ matrices. In particular, when
$$
A_1=\pmatrix{1\\ &1&1\\ &0&-1}\quad\text{and}\quad B_1=\pmatrix{1\\ &0&1\\ &1&0},
$$
we have $A_1^2=B_1^2=I_3$ but $A_1B_1\ne B_1A_1$ over any field.
If you think about how to prove the statement in question, you will see what is special about $2\times2$ matrices. By Cayley-Hamilton theorem, we have $A^2-\operatorname{tr}(A)A+\det(A)I_2=0$ and the analogous holds for $B$. Therefore, the condition $A^2=B^2$ is equivalent to $\operatorname{tr}(A)A-\det(A)I_2=\operatorname{tr}(B)B-\det(B)I_2$. This, when coupled with the further condition $\operatorname{tr}(A)\ne0$ (we actually don’t need both $A$ and $B$ to have zero trace), simply means that $A=pB+qI$ for some scalars $p$ and $q$. Therefore $A$ can be expressed as a polynomial in $B$ and it commutes with $B$.
Does it have a generalisation?
I don’t know any nice ones. We have seen that behind the scene, the given conditions actually effect a direct manipulation of the coefficients of the characteristic polynomial, so as to make $A^2$ and $B^2$ equal to $a_1A+a_0I$ and $B^2=b_1B+b_0I$ respectively for some scalars $a_0,a_1,b_0,b_1$ with $a_1,b_1\ne0$. If you want to generalise this for $n\times n$ matrices, you need to zero out the coefficients of $x^{n-1},x^{n-2},\ldots,x^2$ in the characteristic polynomials of $A$ and $B$ and require that the coefficients of $x$ are nonzero. These conditions can be expressed in terms of the traces of the powers of $A$ and $B$ (thanks to Newton’s identities), but it is hard to conceal the solution trick.
If you don’t insist on generalising the exact problem statement at hand, there is a problem in a similar style, whose result holds in every dimension:
Let $A,B,C$ be three real or complex $n\times n$ matrices such that $AC=CA$ and $BC=CB$. If $(aA+bB+C)^n=AB-BA$ for some pair of scalars $(a,b)\ne(0,0)$, then $AB=BA$.
(See here for a proof of a similar statement.) Yet I don’t think any average high school student is able to solve this.
Best Answer
Conceptually, consider the linear transformations defined by rotating the space through an angle $\theta$, flipping one of the coordinate axes, and rotating back. This is clearly an involution, and since it is a linear transformation, there is a matrix representing it. They are also all distinct, because for each $\theta$, there is a different unit vector that gets flipped around exactly (it's the one that gets rotated to the flipped axis in the first step). Since there are infinitely many angles $\theta$, there are infinitely many such transforms.
For a more general version of this argument, let $D$ be a diagonal matrix whose nonzero elements are all $\pm 1$ that is not $I$ or $-I$ and let $O$ be an orthogonal matrix. Then $$ T = O^TDO $$ is an involution, and this involution will in general be distinct for distinct $O$.