  We present algorithms to solve coupled systems of linear differential
equations, arising in the calculation of massive Feynman diagrams with local
operator insertions at 3-loop order, which do {\it not} request special choices
of bases. Here we assume that the desired solution has a power series
representation and we seek for the coefficients in closed form. In particular,
if the coefficients depend on a small parameter $\ep$ (the dimensional
parameter), we assume that the coefficients themselves can be expanded in
formal Laurent series w.r.t.\ $\ep$ and we try to compute the first terms in
closed form. More precisely, we have a decision algorithm which solves the
following problem: if the terms can be represented by an indefinite nested
hypergeometric sum expression (covering as special cases the harmonic sums,
cyclotomic sums, generalized harmonic sums or nested binomial sums), then we
can calculate them. If the algorithm fails, we obtain a proof that the terms
cannot be represented by the class of indefinite nested hypergeometric sum
expressions. Internally, this problem is reduced by holonomic closure
properties to solving a coupled system of linear difference equations. The
underlying method in this setting relies on decoupling algorithms, difference
ring algorithms and recurrence solving. We demonstrate by a concrete example
how this algorithm can be applied with the new Mathematica package
\texttt{SolveCoupledSystem} which is based on the packages \texttt{Sigma},
\texttt{HarmonicSums} and \texttt{OreSys}. In all applications the
representation in $x$-space is obtained as an iterated integral representation
over general alphabets, generalizing Poincar\'{e} iterated integrals.
