  We present a numerical implementation of the recently developed
unconditionally convergent representation of general Heun functions as integral
series. We produce two codes in Python available for download, one of which is
especially aimed at reproducing the output of Mathematica's HeunG function. We
show that the present code compares favorably with Mathematica's HeunG and with
an Octave/Matlab code of Motygin, in particular when the Heun function is to be
evaluated at a large number of points if less accuracy is sufficient. We
suggest further improvements concerning the accuracy and discuss the issue of
singularities.
