A magic square matrix $M$ is a square matrix with real entries such that the sum of the entries in each column, each row, and each main diagonal is the same.
The problem is to characterize all $3 \times 3$ magic square matrices $M$ such that $M^2$ is also a magic square matrix.
Supposedly, there is an elegant way to do this besides writing out variables and bashing out the multiplication, but I haven't yet found such a solution.
It can be shown that a generic $3\times 3$ magic matrix is of the form: $$\left[\begin{array}{ccc}a&b&c\\d&e&f\\g&h&i\end{array}\right]=\left[\begin{array}{ccc}e+h-c&2e-h&c\\2c-h&e&2e-2c+h\\2e-c&h&e-h+c\end{array}\right].$$ Therefore, the square of a magic matrix can be expressed as: $$\left[ \begin{array}{ccc} e^2+2 h^2+4 c e-4ch & 4 e^2-h^2-2 c e+2 c h & 4 e^2-h^2-2 c e+2 c h \\ 4 e^2-h^2-2 c e+2 c h & e^2+2 h^2+4 c e-4ch & 4 e^2-h^2-2 c e+2 c h \\ 4 e^2-h^2-2 c e+2 c h & 4 e^2-h^2-2 c e+2 c h & e^2+2 h^2+4 c e-4ch \\ \end{array} \right].$$ We want this to be of the upper form, which gives only one equation: $$(2 c - e - h) (e - h)=0.$$ Thus, $h=e$ or $h=2c-e$ are our solutions. The requested matrices are: $$\left[\begin{array}{ccc}2e-c&e&c\\2c-e&e&3e-2c\\2e-c&e&c\end{array}\right],\left[\begin{array}{ccc}c&3e-2c&c\\e&e&e\\2e-c&2c-e&2e-c\end{array}\right],$$ for any $c,e\in\mathbb{R}$.