We propose a new model of the TeV to PeV Galactic diffuse gamma-ray emission, that takes into account the disreteness of cosmic-ray injection sites, by implementing numerically the analytical solutions of the diffusion equation obtained by means of the Green’s functions formalism. As an application of our model, we compute the VHE diffuse gamma-ray emission for two different cosmic-ray diffusion processes and demonstrate that our approach leads to a more realistic description of the gamma-ray sky. In particular, our simulations show that for some relevant regions of the parameters space, such as considering that only a small subset of the cosmic-ray sources are PeVatrons, the discretness of the cosmic-ray sources can no longer be ignored and induces some clumps on the gamma-ray sky maps. In the light of these new results, we expect that near-future observations will provide relevant constraints on the cosmic-ray production and
propagation mechanisms.