In this study, a collocation method based on the Bernoulli polynomials is presented to nd approximate solutions of a system of high-order linear Volterra Integral Equations (VIEs) with variable coecients. In fact, the approximate solution of the problem in the truncated Bernoulli series form is obtained by this method. In addition, the method is presented with error and stability analysis. To show the accuracy and eciency of the method, numerical examples are implemented and the comparisons are given by the other methods.