Thanks to Higham I know that $A f(BA) = f(AB) A$ for any two matrices whose sizes are compatible.

Now I believe that $A (BA)^D = (AB)^D A$, even though the Drazin inverse is not the same function (polynomial?) for $AB$ as for $BA$.

I have validated this relationship via numerical experiments with random matrices, I just can't $prove$ it.

Can you prove (or disprove) it?