In this study, Taylor matrix algorithm is designed for the approximate solution of linear and non-linear differential equation systems. The algorithm is essentially based on the expansion of the functions in differential equation systems to Taylor series and substituting the matrix forms of these expansions into the given equation systems. Using the Mathematica program, the matrix equations are solved and the unknown Taylor coefficients are found approximately. The presented numerical approach is discussed on samples from various linear and non-linear differential equation systems as well as stiff systems. The computational data are then compared with those of some earlier numerical or exact results. As a result, this comparison demonstrates that the proposed method is accurate and reliable.