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$.
 A: Let $w$ be a primitive cubic root of unity and $X=A+wB$. If $A^2+B^2=AB$, then
$$
X\overline{X}=(A+wB)(A+w^2B)=wBA+(1+w^2)AB=w(BA-AB).
$$
Therefore $w^2\det(BA-AB)=|\det(X)|^2$ is real and $\det(BA-AB)$ must be zero. It follows that $X$ is singular. As it is also $2\times2$, we may write $X=uv^\ast$ for some $u,v\in\mathbb C^2$. Hence
$$
|v^\top u|^2=\operatorname{tr}\left(u(\overline{v^\top u})v^\top\right)
=\operatorname{tr}(X\overline{X})=w\operatorname{tr}(AB-BA)=0
$$
and $v^\top u=0$. Consequently, $w(AB-BA) = X\overline{X} = u(\overline{v^\top u})v^\top = 0$, i.e. $AB=BA$.
A: $\textbf{Proposition 1}$. Let $A,B\in M_2(\mathbb{C})$ s.t. $A^2+B^2=AB$. Then $A,B$ have a common eigenvector.
$\textbf{Proof}$. Let $A=X+Y,B=-jX-j^2Y$ where $j=\exp(2i\pi/3)$. Then 
$XY=j^2YX, AB-BA=(j-j^2)(XY-YX)$.
Note that $\det(XY)=j\det(YX)$ implies that $\det(XY)=0$. Clearly $\det(AB-BA)=0$ and $rank(AB-BA)\leq 1$. That implies that $A,B$ have a common eigenvector. $\square$
$\textbf{Proposition 2}$. Let $A,B\in M_2(\mathbb{R})$ s.t. $A^2+B^2=AB$. Then $AB=BA$.
$\textbf{Proof}$. Let $u\in\mathbb{C}^2\setminus\{0\}$ s.t. $Au=\lambda u,Bu=\mu u$. 


*

*If we cannot choose $u$ in $\mathbb{R}^2\setminus\{0\}$ then necessarily $\lambda,\mu\notin\mathbb{R}$ and $A\overline{u}=\overline{\lambda}\overline{u},B\overline{u}=\overline{\mu}\overline{u}$. Since $\{u,\overline{u}\}$ is a basis of $\mathbb{C}^2$, we deduce that $AB=BA$.

*If we can choose $u$ in $\mathbb{R}^2\setminus\{0\}$ then necessarily $\lambda,\mu\in\mathbb{R}$ and the eigenvalues of $A,B$ are all real. Thus, we may assume that $A,B$ are upper-triangular real matrices in the form 
$A=\begin{pmatrix} a&b\\0&c\end{pmatrix},B=\begin{pmatrix}d&e\\0&f\end{pmatrix}$ where $f\in\{-jc,-j^2c\},d\in\{-ja,-j^2a\},b.tr(A)+e.tr(B)=ae+bf$. 
We conclude that $a=c=d=f=0$ and $AB=BA$. $\square$
Note that there are solutions over $\mathbb{C}$ that do not commute.
A: 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$.
