We develop and apply the product integration method to a large class of linear weakly singular Volterra
systems. We show that under certain sufficient conditions this
method converges. Numerical implementation of the method is
illustrated by a benchmark problem originated from heat
conduction.