In this paper, we propose a class of methods to solve the parabolic Volterra integro-differential equations with bounded and unbounded domains. More precisely, we change the parabolic Volterra integro-differential equations to well-posed linear and nonlinear dynamical systems. Then, the obtained systems are solved by using a new class of algorithms consisting linear multi-step formulas in which these schemes are constructed through the hybrid of Gergory's formula, finite difference and multi-step methods. Error bounds are derived in both bounded and unbounded domains. Some numerical examples are then presented to illustrate the efficiency and accuracy of the proposed methods. Furthermore, stability and convergence of proposed methods are established and we denote the numerical simulations. Moreover, some tests are conducted on data with measurement noise to consider the performance of the proposed methods.