Solving differential equations for Feynman integrals by expansions near singular points

We describe a strategy to solve differential equations for Feynman integrals by powers series expansions near singular points and to obtain high precision results for the corresponding master integrals. We consider Feynman integrals with two scales, i.e. nontrivially depending on one variable. The corresponding algorithm is oriented at situations where canonical form of the differential equations is impossible. We provide a computer implementation of our algorithm in a simple example of four-loop generalized sun-set integrals with three equal non-zero masses. Our code provides values of the master integrals at any given point on the real axis with a required accuracy and a given order of expansion in the regularization parameter $\epsilon$.