  This article is devoted to Feller's diffusion equation which arises naturally
in probabilities and physics (e.g. wave turbulence theory). If discretized
naively, this equation may represent serious numerical difficulties since the
diffusion coefficient is practically unbounded and most of its solutions are
weakly divergent at the origin. In order to overcome these difficulties we
reformulate this equation using some ideas from the Lagrangian fluid mechanics.
This allows us to obtain a numerical scheme with a rather generous stability
condition. Finally, the algorithm admits an elegant implementation and the
corresponding Matlab code is provided with this article under an open source
license.
