$A^2+B^2=AB$ and $BA-AB$ is non-singular

The answer is yes. We shall construct a concrete example below. First, suppose $A$ and $B$ are invertible and let $C=AB^{-1}$. From $A^2+B^2=AB$ we get $C + A^{-1}C^{-1}A = I$ and $A^{-1}CA=(I-C)^{-1}$. Therefore the spectrum of $C$ is invariant under the transformation $z\mapsto \frac1{1-z}$.

Let $w$ be the cubic root of unity. If $A^2+B^2=AB$, then $(A+wB)(A+w^2B)=w(BA-AB)$. So, if we want $BA-AB$ to be invertible, neither $-w$ nor $-w^2$ can be an eigenvalue of $C$. As $-w$ and $-w^2$ are fixed points of $z\mapsto\frac1{1-z}$, the spectrum of $C$ has to be splittable into three cycles of the form $\{z,\ \frac1{1-z},\ \ 1-\frac1z\}$, where $z\ne-w,-w^2$, and for every feasible pair of invertible matrices $A,B$, their dimensions must be divisible by 3.

Now we can construct an example on $M_3(\mathbb C)$. Take $z=2$ and $C=\operatorname{diag}(2, -1, \frac12)$. So $(I-C)^{-1}=\operatorname{diag}(-1, \frac12, 2)$. To make $A^{-1}CA=(I-C)^{-1}$ or $A^2+B^2=AB$, we pick $A=\pmatrix{0&0&1\\ 1&0&0\\ 0&1&0}$ and $B=C^{-1}A=\pmatrix{0&0&\frac12\\ -1&0&0\\ 0&2&0}$. One can readily verify that $BA-AB=\pmatrix{0&-\frac32&0\\ 0&0&-\frac32\\ 3&0&0}$ is invertible.

We consider the equations
(1) $A^2+B^2=AB$,
(2) $A^2+B^2=2AB$.
A couple $(A,B)$, solution of (2), is simultaneously triangularizable, and, moreover, when $n=2$, $AB=BA$.

Eq. (1) has not the second property: indeed $$A=\begin{pmatrix}0&-1\\1&-1\end{pmatrix},\ B=\begin{pmatrix}\dfrac{1-i\sqrt{3}}{2}&0\\\dfrac{3+i\sqrt{3}}{2}&-1\end{pmatrix}$$ is a solution s.t. $AB\not= BA$. Note that $\mathrm{spectrum}(A)=\{e^{\pm 2i\pi/3}\}$ and $A^3=I_2$, $B^3=-I_2$; moreover $(AB-BA)^2=0$. I think that Eq. (1) has the first property, or at least, $AB-BA$ is nilpotent, but I don't know how to prove it. We see, in the following, that, a solution, in general position, satisfies $AB=BA$.

Proposition. We consider the equation in the unknown $X$: $X^2-AX+A^2=0$. Let $\mathrm{spectrum}(A)=(\lambda_i)_i$ denote the complete spectrum of $A$. If, for every $i\not= j$ $\lambda_i^3\not=\lambda_j^3$, then any solution $X$ satisfies $AX=XA$.

Proof. Since $A$ is diagonalizable we may assume $A=\mathrm{diag}(\lambda_i)$. Thus $A^3=\mathrm{diag}(\lambda_i^3)$ has distinct eigenvalues; since $X$ and $A^3$ commute, $X$ is a diagonal matrix and we are done.

EDIT. There are solutions of $A^3+B^3=0$ s.t. $AB-BA$ is invertible; for example, choose $A=\begin{pmatrix}0&1&0\\0&0&1\\0&0&0\end{pmatrix},\ B=\begin{pmatrix}0&0&0\\1&0&0\\0&-1&0\end{pmatrix}$; of course $A^2+B^2\not= AB$.