This study mainly investigates new techniques for obtaining numerical solutions of time-fractional diffusion equations. The fractional derivative term is represented in the Lagrange operational sense. First, we describe the temporal direction of the considered model using the Legendre orthogonal polynomials. Moreover, to archive a full discretization approach a type of nonlocal method has been applied that is known as the nonlocal peridynamic differential operator (PDDO). The PDDO is based on the concept of peridynamic (PD) interactions by proposing the PD functions orthogonal to each term in the Taylor Series Expansion (TSE) of a field variable. The PDDO for numerical integration uses the vicinity of each point (referred to as the horizon, which does not need background mesh). The PDDO is exclusively described in terms of integration (summation) throughout the interaction domain. As a result, it is unsusceptible to singularities caused by discontinuities. It does, however, need the creation of PD functions at each node. We numerically investigate the stability, the convergence of the scheme, which verifies the validity of the proposed method. Numerical results show the simplicity and accuracy of the presented method.