We address the influences of heterogeneity and complex geology of the matrix of fractured rocks on transient flow. Fractional constitutive flux laws reflect the stochastic framework that we consider. We model transient diffusion in both the matrix and fracture systems in terms of a continuous time random walk. Our procedure is particularly suited to address a complex geology that may exist on a number of scales and which may include dead ends and discontinuities. The performance of a horizontal well produced through arbitrarily located, multiple hydraulic fractures with distinct properties (length, width, permeability) forms the basis for our thesis. The pressure distribution in a rectangular drainage region where the well may be placed arbitrarily is expressed in terms of the Laplace transformation. The required solutions are obtained numerically. The focus of our work is on long-term behaviors for production at a constant pressure. In addition to numerical solutions, asymptotic results that provide information on the structure of the solutions are presented. Agreement between the asymptotic and numerical solutions is excellent. We show that long-term responses are governed by two distinct, two-parameter Mittag-Leffler functions and are an outcome of the complexities we desire to model in both the matrix and fracture systems. As a consequence, power-law behaviors that reflect the heterogeneity inherent in the system define long-time expectations. We show that the new solutions we derive do reduce to those of classical diffusion; that is, results corresponding to classical diffusion are a subset of the new results obtained here. Our results are particularly suited to model transient flow in shale reservoirs.