If your matrices have exact entries (for example in ℚ or a finite field), you can consider the polynomial ring $R = K[x_{i,j} | 0≤i,j\lt n]$ in $n×n$ variables, i.e. the variables $x_{i,j}$ represent entries of an arbitrary matrix $X$.
sage: K = QQ
sage: n = 3
sage: S = [matrix(K, [[0,2,0],[0,0,-1],[1,0,0]])]
sage: R = PolynomialRing(K, ['x_%d_%d' % (a, b) for a in range(n) for b in range(n)])
sage: X = matrix(R, n, R.gens()); X
[x_0_0 x_0_1 x_0_2]
[x_1_0 x_1_1 x_1_2]
[x_2_0 x_2_1 x_2_2]
The condition $X Y - Y X = 0$ (entrywise) imposes linear relations on the variables $x_{i,j}$. These form an ideal $J$.
sage: J = R.ideal(flatten([(X*Y-Y*X).list() for Y in S]))
sage: J.dimension()
3
As the ideal is generated by linear forms, this (Krull) dimension of the ideal should agree with the vector space dimension you are after (I think). It should agree with the length of
sage: J.normal_basis(1)
[x_2_2, x_2_1, x_2_0]
which gives the variables that parametrize $X$ modulo $J$. We can map $X$ to the quotient ring:
sage: Q = R.quotient(J)
sage: X.change_ring(Q)
[ x_2_2bar 2*x_2_0bar -x_2_1bar]
[-1/2*x_2_1bar x_2_2bar -x_2_0bar]
[ x_2_0bar x_2_1bar x_2_2bar]
The result indeed only depends on these three parameters. From there, you could construct e.g. a basis by substituting these three variables by elements from $K$.
Alternatively, to avoid $Q$ altogether, you can compute $X$ modulo $J$ like this:
sage: X.apply_map(J.reduce)
[ x_2_2 2*x_2_0 -x_2_1]
[-1/2*x_2_1 x_2_2 -x_2_0]
[ x_2_0 x_2_1 x_2_2]
