This paper contributes a new matrix method for solving systems of high-order linear differential–difference equations with variable coefficients under given initial conditions. On the basis
of the presented approach, the matrix forms of the Euler polynomials and their derivatives are constructed, and then by substituting the collocation points into the matrix forms, the fundamental
matrix equation is formed. This matrix equation corresponds to a system of linear algebraic equations. By solving this system, the unknown Euler coefficients are determined. Some illustrative
examples with comparisons are given. The results demonstrate reliability and efficiency of the proposed method.