Positive: The routine computes the function in regions where the Taylor series fails.
Negative: Computations are relatively slow and not completely reliable, for example, alpha=1/2, beta=1 and z near the imaginary axis (probably directly on the imaginary axis) yield completely wrong results (see the phase portrait at http://www.mathe.tu-freiberg.de/~wegert/MittagLefflerBug/Mittag4.png)
I have one more remark. I found this article on internet:
and it seemed to me that your code is based on it. Probably a reference in your code to this paper (or the theory behind your code), can help other people understand your code.