The non-steady transport and decay of pollutants in drinking water distribution netwolks is traditionally modeled with pure advection models. This paper presents an eulerian-lagrangian numerical solution for the advection-diffusion equatior in wate supply networks Because of the diffusion term, the numerical scheme produces a large linear system of equations. A technique that uses numerically computed Green's functions and nodal mass balance considerations is proposed to decompose this large system in three tri-diagona systems for each pipe and one (much smaller) system for the concentration at the pipe junctions In this way the system is efficiently solved and the model can be applied to large networks at reasonable computer cost. The model is applied to simulate the fluoride propagation in a real network for which published data of field measurements and simulation with the U.S.E.P.A. EPANET model are available. Comparisons are presented in the paper that clearly show that while in pipes with high and medium values of flow velocity the two models give very similar results, in the pipes with low values of flow velocity the new advection-diffusion mode more closely predicts the concentration evolution, provided an appropriate value for the dispersion coefficient is used.