  The computation of the Mittag-Leffler (ML) function with matrix arguments,
and some applications in fractional calculus, are discussed. In general the
evaluation of a scalar function in matrix arguments may require the computation
of derivatives of possible high order depending on the matrix spectrum.
Regarding the ML function, the numerical computation of its derivatives of
arbitrary order is a completely unexplored topic; in this paper we address this
issue and three different methods are tailored and investigated. The methods
are combined together with an original derivatives balancing technique in order
to devise an algorithm capable of providing high accuracy. The conditioning of
the evaluation of matrix ML functions is also studied. The numerical
experiments presented in the paper show that the proposed algorithm provides
high accuracy, very often close to the machine precision.
