Heavy-ion collisions at center-of-mass energies between 1 and 100 GeV/nucleon are essential to understand the phase diagram of QCD and search for its critical point. At these energies the net baryon density of the system can be high, and simulating its evolution becomes an indispensable part of theoretical modeling. We here present the (3+1)-dimensional diffusive relativistic hydrodynamic code BEShydro which solves the equations of motion of second-order Denicol–Niemi–Molnar–Rischke (DNMR) theory, including bulk and shear viscous currents and baryon diffusion currents. BEShydro features a modular structure that allows to easily turn on and off baryon evolution and different dissipative effects and thus to study their physical effects on the dynamical evolution individually. An extensive set of test protocols for the code, including several novel tests of the precision of baryon transport that can also be used to test other such codes, is documented here and supplied as a permanent part of the code package. Program Title:BEShydro Program Files doi:http://dx.doi.org/10.17632/xwywvb2psm.1 Licensing provisions: GPLv3 Programming language: C++ External routines/libraries: GNU Scientific Library (GSL) Nature of problem: (3+1)-dimensional dynamical evolution of hot and dense matter created in relativistic heavy-ion collisions using second-order dissipative relativistic fluid dynamics, including evolution of net baryon number and its dissipative diffusion current. Solution method: Runge–Kutta Kurganov–Tadmor algorithm.