This paper deals with the computation of the two parameter Mittag-Leffler function of operators by exploiting its Stieltjes integral representation and then by using a single exponential transform together with the sinc rule. Whenever the parameters of the function do not allow this representation, we resort to the Dunford-Taylor one. The error analysis is kept in the framework of unbounded accretive operators in order to make it a useful tool for the solution of fractional differential equations. The theory is also used to design a rational Krylov method.