We propose algorithms for the computation of the first N terms of a vector (or a full basis) of power series solutions of a linear system of differential equations at an ordinary point, using a number of arithmetic operations that is quasi-linear with respect to N. Similar results are also given in the nonlinear case. This extends previous results obtained by Brent and Kung for scalar differential equations of order 1 and 2.