The main aim of this paper is to construct an efficient highly accurate numerical scheme to solve a class of one and two-dimensional parabolic integro-fractional differential equations. The high order $L2$-$1_\sigma$ scheme is taken into account to discretize the time-fractional operator on a uniform mesh. To solve the one-dimensional problem, second-order discretizations are used to approximate the spatial derivatives whereas, a repeated quadrature rule based on trapezoidal approximation is employed to discretize the integral operator. In the case of two-dimensional problem, first, we make the semi-discretization of the proposed model based on the $L2$-$1_\sigma$ scheme for the fractional operator, and composite trapezoidal approximation for the integral part. Then, the spatial derivatives are approximated by the two-dimensional Haar wavelet. The stability and convergence analysis is carried out for both models. The experimental evidence proves the strong reliability of the present methods. Further, the obtained results are compared with some existing methods through several graphs and tables, and it is shown that the proposed methods not only have better accuracy but also produce less error in comparison with the $L1$ scheme.