In this paper, a numerical matrix method based on collocation points is presented for the approximate solution of the systems of high-order linear Fredholm integro-differential equations with variable coefficients under the mixed conditions in terms of the Bessel polynomials. Numerical examples are included to demonstrate the validity and the applicability of the technique and also the results are compared with the different methods. The results show the efficiently and the accuracy of the present work. All of the numerical computations have been performed on a PC using some programs written in MATLAB v7.6.0 (R2008a).