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.