A search for heavy neutral leptons (HNLs) decaying in the CMS muon system is presented. A data sample is used corresponding to an integrated luminosity of 138  fb⁻¹ of proton-proton collisions at √s=13  TeV, recorded at the CERN LHC in 2016–2018. Decay products of long-lived HNLs could interact with the shielding materials in the CMS muon system and create hadronic and electromagnetic showers detected in the muon chambers. This distinctive signature provides a unique handle to search for HNLs with masses below 4 GeV and proper decay lengths of the order of meters. The signature is sensitive to HNL couplings to all three generations of leptons. Candidate events are required to contain a prompt electron or muon originating from a vertex on the beam axis and a displaced shower in the muon chambers. No significant deviations from the standard model background expectation are observed. In the electron (muon) channel, the most stringent limits to date are set for HNLs in the mass range of 2.1–3.0 (1.9–3.3) GeV, reaching mixing matrix element squared values as low as 8.6⁢(4.6)×10⁻⁶.