An algorithm to find a solution
of differential equations for master integrals
in the form of an $\epsilon$-expansion series with numerical coefficients is applied
to 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 20000 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.