We consider systems of higher-order linear differential equations having a singularity at the origin. We develop an efficient method that computes their regular formal solutions. Our algorithm handles the system directly without transforming it into a system of first order but higher size. We generalize the approach proposed by Poole in the scalar case consisting in arranging the solution as a Frobenius series in the variable $x$ which coefficients are polynomials in $\ln(x)$. The problem is then reduced to computing polynomial solutions of non-homogeneous linear differential systems with constant matrix coefficients. This latter problem is solved using the properties of matrix polynomials, the determinant of which plays the same role as indicial polynomial in the scalar case. We give the arithmetic complexity of our algorithm and compare our implementation in Maple to existing methods. This is a joint work with M. Barkatou and T.Cluzeau.