The paper introduces an acurrate numerical appraoch based on orthonormal Bernoulli polynomials for solving parabolic partial integro- differential equations (PIDEs). This type of equations arises in physics and engineering. Some operational matrix are given for these polynomials and are also used to obtain the numerical solution. By this approach, the problem is transformed into a nonlinear algebraic system. Convergence analysis is given and some experiment tests are studied to examine the good accuracy of the numerical algorithm, the proposed technique is compared with some other well known methods.