The problem of mass diffusion in layered systems has relevance to applications in different scientific disciplines, e.g., chemistry, material science, soil science, and biomedical engineering. The mathematical challenge in these type of model systems is to match the solutions of the time-dependent diffusion equation in each layer, such that the boundary conditions at the interfaces between them are satisfied. As the number of layers increases, the solutions may become increasingly complicated. Here, we describe an alternative computational approach to multi-layer diffusion problems, which is based on the description of the overdamped Brownian motion of particles via the underdamped Langevin equation. In this approach, the probability distribution function is computed from the statistics of an ensemble of independent single particle trajectories. To allow for simulations of Langevin dynamics in layered systems, the numerical integrator must be supplemented with algorithms for the transitions across the discontinuous interfaces. Algorithms for three common types of discontinuities are presented: (i) A discontinuity in the friction coefficient, (ii) a semi-permeable membrane, and (iii) a step-function chemical potential. The general case of an interface where all three discontinuities are present (Kedem-Katchalsky boundary) is also discussed. We demonstrate the validity and accuracy of the derived algorithms by considering a simple two-layer model system and comparing the Langevin dynamics statistics with analytical solutions and alternative computational results.