Diffusion processes intended to model the continuous state space limit of birth–death processes, chemical reactions, and other discrete particle systems often involve multiplicative noise where the diffusion vanishes near one (or more) of the state space boundaries. Standard direct numerical simulation schemes for the associated stochastic differential equations run the risk of “overshooting”, i.e., of varying outside the meaningful state space domain where simple analytic expressions for the diffusion coefficient may take on unphysical (negative or complex) values. We propose a simple scheme to overcome this problem and apply it to an exactly soluble stochastic ordinary differential equation (SODE), and to a related parabolic stochastic partial differential equation (SPDE) that admits exact analytic solution for the stationary correlation function. Armed with these analytic benchmark solutions, we demonstrate that the scheme produces approximate solutions for the SODE with distributions that display first-order convergence in the Wasserstein metric. For the SPDE, the scheme produces first order convergence for the stationary correlation function in L2.