  We outline a new algorithm to solve coupled systems of differential equations
in one continuous variable $x$ (resp. coupled difference equations in one
discrete variable $N$) depending on a small parameter $\epsilon$: given such a
system and given sufficiently many initial values, we can determine the first
coefficients of the Laurent-series solutions in $\epsilon$ if they are
expressible in terms of indefinite nested sums and products. This systematic
approach is based on symbolic summation algorithms in the context of difference
rings/fields and uncoupling algorithms. The proposed method gives rise to new
interesting applications in connection with integration by parts (IBP) methods.
As an illustrative example, we will demonstrate how one can calculate the
$\epsilon$-expansion of a ladder graph with 6 massive fermion lines.
