We propose a nodal high-order discontinuous Galerkin (DG) method for coupled wave propagation in heterogeneous elastoplastic soil columns. We solve the elastodynamic system written in velocity-strain formulation considering simultaneously the three components of motion. A 3D Iwan elastoplastic model for dry soils under dynamic loading is introduced in the DG method (DG-3C) based on centered fluxes and a low storage Runge–Kutta time scheme. Unlike the linear case, the nonlinear material behavior results in coupling effects between different components of motion. We focus on the choice of the numerical flux for nonlinear heterogeneous media and the classical centered flux is well adapted in terms of stability and direct implementation in 3D. Several numerical applications are studied using a soil column with realistic mechanical properties, and considering successively synthetic and borehole seismic recordings recorded at the Volvi test site (northern Greece).