We elaborate on the method of differential equations for evaluating Feynman integrals. We focus on systems of equations for master integrals having a linear dependence on the dimensional parameter. For these systems we identify the criteria to bring them in a canonical form, recently identified by Henn, where the dependence of the dimensional parameter is disentangled from the kinematics. The determination of the transformation and the computation of the solution are obtained by using Magnus and Dyson series expansion. We apply the method to planar and non-planar two-loop QED vertex diagrams for massive fermions, and to non-planar two-loop integrals contributing to 2 → 2 scattering of massless particles. The extension to systems which are polynomial in the dimensional parameter is discussed as well.