In this study, a collocation method based on the Bernoulli polynomials is presented to find approximate solutions of a system of high-order linear Volterra integral equations (VIEs) with variable coefficients. 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 the efficiency of the method, numerical examples are implemented and the comparisons are given by the other methods.