Evaluating `elliptic' master integrals at special kinematic values: using differential equations and their solutions via expansions near singular points

This is a sequel of our previous paper where we described an algorithm to find a solution of differential equations for master integrals in the form of an $\epsilon$-expansion series with numerical coefficients. The algorithm is based on using generalized power series expansions near singular points of the differential system, solving difference equations for the corresponding coefficients in these expansions and using matching to connect series expansions at two neighboring points. Here we use our algorithm and the corresponding code for our example of four-loop generalized sunset diagrams with three massive and two massless propagators, in order to obtain new analytical results. We analytically evaluate the master integrals at threshold, $p^2=9 m^2$, in an expansion in $\epsilon$ up to $\epsilon^1$. With the help of our code, we obtain numerical results for the threshold master integrals in an $\epsilon$-expansion with the accuracy of 6000 digits and then use the PSLQ algorithm to arrive at analytical values. Our basis of constants is build from bases of multiple polylogarithm values at sixth roots of unity.