  We present an algorithm which allows to solve analytically linear systems of
differential equations which factorize to first order. The solution is given in
terms of iterated integrals over an alphabet where its structure is implied by
the coefficient matrix of the differential equations. These systems appear in a
large variety of higher order calculations in perturbative Quantum Field
Theories. We apply this method to calculate the master integrals of the
three--loop massive form factors for different currents, as an illustration,
and present the results for the vector form factors in detail. Here the
solution space emerging is given by the cyclotomic harmonic polylogarithms and
their associated special constants. No special basis representation of the
master integrals is needed. The algorithm can be applied as well to more
general cases factorizing at first order, which are based on more general
alphabets, iterated integrals and associated constants.
