Departing from a general stochastic model for a moving boundary problem, we consider the density function of probability for the first passing time. It is well known that the distribution of this random variable satisfies a problem ruled by an advection–diffusion system for which very few solutions are known in exact form. The model considers also a deterministic source, and the coefficients of this equation are functions with sufficient regularity. A numerical scheme is designed to estimate the solutions of the initial-boundary-value problem. We prove rigorously that the numerical model is capable of preserving the main characteristics of the solutions of the stochastic model, that is, positivity, boundedness and monotonicity. The scheme has spatial symmetry, and it is theoretically analyzed for consistency, stability and convergence. Some numerical simulations are carried out in this work to assess the capability of the discrete model to preserve the main structural features of the solutions of the model. Moreover, a numerical study confirms the efficiency of the scheme, in agreement with the mathematical results obtained in this work.