Reaction-diffusion systems are mathematical models that describe how the concentrations of substances distributed in space change under the influence of local chemical reactions, and diffusion which causes the substances to spread out in space. The classical represen- tation of a reaction-diffusion system is given by semi-linear parabolic partial differential equations, whose solution predicts how diffusion causes the concentration field to change with time. This change is proportional to the diffusion coefficient. If the solute moves in a homogeneous system in thermal equilibrium, the diffusion coefficients are constants that do not depend on the local concentration of solvent and solute. However, in non- homogeneous and structured media the assumption of constant intracellular diffusion co- efficient is not necessarily valid, and, consequently, the diffusion coefficient is a function of the local concentration of solvent and solutes. In this paper we propose a stochastic model of reaction-diffusion systems, in which the diffusion coefficients are function of the local concentration, viscosity and frictional forces. We then describe the software tool Redi (REaction-DIffusion simulator) which we have developed in order to implement this model into a Gillespie-like stochastic simulation algorithm. Finally, we show the ability of our model implemented in the Redi tool to reproduce the observed gradient of the bicoid protein in the Drosophila Melanogaster embryo. With Redi, we were able to simulate with an accuracy of 1% the experimental spatio-temporal dynamics of the bicoid protein, as recorded in time-lapse experiments obtained by direct measurements of transgenic bicoid- enhanced green fluorescent protein.