In this paper, a numerical matrix method, which is based on collocation points, is presented for the approximate solution of a system of high-order linear differential-difference equations with variable coefficients under the mixed conditions in terms of Bessel polynomials. Numerical examples are included to demonstrate the validity and applicability of the technique and comparisons are made with existing results. The results show the efficiency and accuracy of the present work. All of the numerical computations have been performed on computer using a program written in MATLAB v7.6.0 (R2008a). © 2011 Verlag der Zeitschrift für Naturforschung, Tübingen.