Consider a matrix
$$A = \begin{pmatrix} a & b & c \\ d & e & f \\ g & h & i \end{pmatrix}.$$
If $A$ is a magic square, then the sum of the rows must be equal, hence
\begin{align*}
a + b + c &= d + e + f \\
a + b + c &= g + h + i.
\end{align*}
This sum must also be the sum of the columns:
\begin{align*}
a + b + c &= a + d + g \\
a + b + c &= b + e + h \\
a + b + c &= c + f + i.
\end{align*}
Finally, it must be the sum of the two diagonals:
\begin{align*}
a + b + c &= a + e + i \\
a + b + c &= c + e + g.
\end{align*}
This gives us $8$ homogeneous linear equations, which we can write as a matrix equation like so:
$$\begin{pmatrix}
1 & 1 & 1 & -1 & -1 & -1 & 0 & 0 & 0 \\
1 & 1 & 1 & 0 & 0 & 0 & -1 & -1 & -1 \\
0 & 1 & 1 & -1 & 0 & 0 & -1 & 0 & 0 \\
1 & 0 & 1 & 0 & -1 & 0 & 0 & -1 & 0 \\
1 & 1 & 0 & 0 & 0 & -1 & 0 & 0 & -1 \\
0 & 1 & 1 & 0 & -1 & 0 & 0 & 0 & -1 \\
1 & 1 & 0 & 0 & -1 & 0 & -1 & 0 & 0 \\
\end{pmatrix}
\begin{pmatrix}a \\ b \\ c \\ d \\ e \\ f \\ g \\ h \\ i\end{pmatrix} =
\begin{pmatrix}0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0\end{pmatrix}.$$
Solve this, using Gaussian elimination, for $a, b, c, d, e, f, g, h, i$ in general, making sure to take into account free parameters. If you write out the solution $A$ in terms of these free parameters, you will automatically have $A$ as a linear combination of a set of magic squares, which form a basis for all $3 \times 3$ magic squares (the parameters will form the coefficients of these linear combinations).
Good luck!