Combining arbitrary order global Pad\'e approximation of the Mittag-Leffler function with its addition formula for a significant accuracy boost

The combination of the global Pad\'e approximation of the Mittag-Leffler function with its addition formula for the case $\alpha<1$ yields significantly higher accuracy results for a given arbitrary order $n$. We present a solution in terms of a Mathematica notebook to determine the general structure of the system of linear equations to be solved, followed by an implementation as a {\tt{C++}} program using the {\tt{Eigen}} template library for linear algebra. For a comparison with contour integral solutions we present an implementation as a {\tt{C++}} program using the {\tt{boost}} library's quadrature package employing the Gauss-Kronrod-method.