<s>[INST] If $A,B \in \mathcal{M_2}(\mathbb{R})$ and $A^2+B^2=AB$, does it follow that $A$ and $B$ commute? Let $A,B \in \mathcal{M_2}(\mathbb{R})$ such that $A^2+B^2=AB$. Is it necessary that $AB=BA$?
I could easily show that such matrices have the property that $(AB-BA)^2=O_2$ (this was actually the question I was asked, then I started wondering if it would be true that the matrices actually commute) by considering the matrix $M=A+\epsilon B$ and computing the determinat  of $M \cdot \overline{M}$($\epsilon \in \mathbb{C}\setminus \mathbb{R}$ is a cubic root of unity), but that's all I got. I tried to find some counterexamples, but I have a hard time finding any matrices with that property.
EDIT: To see that $(AB-BA)^2=O_2$ we do the following : by direct computations 
$$|\det M|^2=\det M \det \overline{M}=\epsilon^2 \det(AB-BA)$$ and this is a real number if and only if $\det(AB-BA)=0$. From the Cayley-Hamilton theorem we now get that $(AB-BA)^2=O_2$.
 [/INST] Proposition: Let $A,B \in \mathcal{M_2}(\mathbb{R})$ such that $A^2+B^2=AB$. Then $AB=BA$.
Disclaimer: This proof is absolutely awful.
Proof. Note that for every real number $\lambda\in\Bbb{R}$ we have
$$(\lambda A)^2+(\lambda B)^2=\lambda^2(A^2+B^2)=\lambda^2(AB)=(\lambda A)(\lambda B).$$
If $\det A\neq0$ then after an appropriate change of basis we have
$$A=\begin{pmatrix}1&0\\0&\lambda\end{pmatrix}
\qquad\text{ or }\qquad
A=\begin{pmatrix}1&\lambda\\0&1\end{pmatrix}
\qquad\text{ or }\qquad
A=\begin{pmatrix}\lambda&-1\\1&\lambda\end{pmatrix},$$
for some $\lambda\in\Bbb{R}$. Let $B=\tbinom{a\ b}{c\ d}$. We treat the three cases separately:


*

*Plugging $A=\tbinom{1\ 0}{0\ \lambda}$ and $B$ into the equation yields the system of equations
\begin{eqnarray*}
1+a^2+bc&=&a,\\
b(a+d)&=&b,\\
c(a+d)&=&\lambda c,\\
\lambda^2+d^2+bc&=&\lambda d.
\end{eqnarray*}
If $bc=0$ then $a^2-a+1=0$ which is impossible. Then $\lambda=a+d=1$ and so $A$ is the identity matrix, so certainly $AB=BA$.

*Plugging $A=\tbinom{1\ \lambda}{0\ 1}$ and $B$ into the equation yields the system of equations
\begin{eqnarray*}
1+a^2+bc&=&a+\lambda c,\\
2\lambda+b(a+d)&=&b+\lambda d,\\
c(a+d)&=& c,\\
1+d^2+bc&=&d.
\end{eqnarray*}
If $bc=0$ then $d^2-d+1=0$ which is imposible. It follows that $a+d=1$ and hence that $2\lambda=\lambda d$. If $\lambda=0$ then $A$ is the identity matrix so certainly $AB=BA$. Otherwise $d=2$ and hence $bc=-3$ and $a=-1$, which leads to $\lambda c=0$, a contradiction.

*Plugging $A=\tbinom{\lambda\ -1}{1\ \hphantom{-}\lambda}$ and $B$ into the equation yields the system of equations
\begin{eqnarray*}
\lambda^2-1+a^2+bc&=&\lambda a-c,\\
-2\lambda+b(a+d)&=&\lambda b-d,\\
2\lambda+c(a+d)&=&a+\lambda c,\\
\lambda^2-1+d^2+bc&=&b+\lambda d.
\end{eqnarray*}
Adding the second and third, and subtracting the last from the first, yields
$$(b+c)(a+d)=\lambda(b+c)+(a-d)
\qquad\text{ and }\qquad
a^2-d^2=\lambda(a-d)-(b+c).$$
Isolating $a-d$ from the former and $b+c$ from the latter and plugging it in shows that
$$a-d=(b+c)(a+d-\lambda)=(a-d)(\lambda-a-d)(a+d-\lambda)=-(a+d-\lambda)^2(a-d),$$
which shows that $a=d$. It follows that $b=-c$ and we are left with the system
\begin{eqnarray*}
\lambda^2-1+a^2-b^2&=&\lambda a-b,\\
-2\lambda+2ab&=&\lambda b-a.
\end{eqnarray*}
The latter shows that $a(2b+1)=\lambda(b+2)$, and multiplying the former by $(b+2)^2$ and substituting and cleaning up yields
$$3a^2(b^2-b+1)=(b+2)^2(b^2-b+1),$$
where $b^2-b+1\neq0$ because $b\in\Bbb{R}$. It follows that $b=-2\pm\sqrt{3}a$ and $\lambda=2a\mp\sqrt{3}$ correspondingly. Then
$$A=\begin{pmatrix}
2a\mp\sqrt{3}&-1\\
1&2a\mp\sqrt{3}
\end{pmatrix}
\qquad\text{ and }\qquad
B=\begin{pmatrix}
a&-2\pm\sqrt{3}{a}\\
2\mp\sqrt{3}{a}&a
\end{pmatrix},$$
and a routine check verifies that again $AB=BA$.
Finally, if $\det A=0$ then after an appropriate change of basis we have
$$A=\begin{pmatrix}\lambda &0\\0&0\end{pmatrix}
\qquad\text{ or }\qquad
A=\begin{pmatrix}0&\lambda\\0&0\end{pmatrix},$$
for some $\lambda\in\Bbb{R}$, where the proposition is trivial if $\lambda=0$. If $\lambda\neq0$ a routine check as before shows that for the first form there is no corresponding matrix $B$ satisfying the identity, and for the second form we see that $B$ must be of the form
$$B=\begin{pmatrix}0&\mu\\0&0\end{pmatrix},$$
for some $\mu\in\Bbb{R}$, which shows that $AB=BA$.
</s>