Nowadays, the conventional Euclidean models are mostly used to describe the behavior of fluid flow through porous media. These models assume the homogeneity of the reservoir, and in naturally fractured reservoir, the fractures are distributed uniformly and use the interconnected fractures assumption. However, several cases such as core, log, outcrop data, production behavior of reservoirs, and the dynamic behavior of reservoirs indicate that the reservoirs have a different behavior other than these assumptions in most cases. According to the fractal theory and the concept of fractional derivative, a generalized dif- fusion equation is presented to analyze the transport in fractal reservoirs. Three outer boundary conditions are investigated. Using exact analytical or semi-analytical solutions for generalized diffusion equation with fractional order differential equation and a fractal physical form, under the usual assumptions, requires large amounts of computation time and may produce inaccurate and fake results for some combinations of parameters. Because of fractionality, fractal shape, and therefore the existence of infinite series, large computation times occur, which is sometimes slowly convergent. This paper provides a computationally efficient and accurate method via differential quadrature (DQ) and gener- alized integral quadrature (GIQ) analyses of diffusion equation to overcome these difficul- ties. The presented method would overcome the imperfections in boundary conditions’ implementations of second-order partial differential equation (PDE) encountered in such problems.