In this paper we develop a numerical scheme based on quadratures to approximate solutions of integro-differential equations involving convolution kernels, \(\nu \), of diffusive type. In particular, we assume \(\nu \) is symmetric and exponentially decaying at infinity. We consider problems posed in bounded domains and in \({\mathbb {R}}\). In the case of bounded domains with nonlocal Dirichlet boundary conditions, we show the convergence of the scheme for kernels that have positive tails, but that can take on negative values. When the equations are posed on all of \({\mathbb {R}}\), we show that our scheme converges for nonnegative kernels. Since nonlocal Neumann boundary conditions lead to an equivalent formulation as in the unbounded case, we show that these last results also apply to the Neumann problem.