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 all the 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.