Recent discoveries of gravitational wave (GW) events most likely originating from black hole (BH) + neutron star (NS) mergers reveal the existence of BH+NS binaries. The formation of BH+NS binaries and their merger rates through isolated binary evolution have been investigated extensively with population synthesis simulations. A detailed stellar evolution modeling of the formation of this population, however, is missing from the literature. In this work, we create the first complete 1D model of more than 30 BH+NS progenitor systems, which are calculated self-consistently until collapse of the iron core with infall velocity exceeding 1000 km s−1. Focusing on the progenitors of BH–NS GW sources, we apply the MESA code starting from a post-common-envelope binary with short orbital period (<1 day) consisting of a BH and a zero-age main-sequence helium star that experiences stable mass transfer. The (ultra)stripped supernova explosion is subsequently modeled using a semianalytic method to reveal final remnant masses and momentum kicks. Three example systems (A, B, and C) eventually evolve into BH+NS binaries with component masses of (M
BH, M
NS) = (8.80, 1.53), (8.92, 1.45), and (5.71, 1.34) M
⊙, respectively. These NS masses could be significantly larger depending on the exact mass cut during the supernova explosion. These BH+NS systems are likely to merge and produce GW events within a Hubble time. System C is a potential progenitor of a GW200115-like event, while Systems A and B are possible candidates for a GW200105-like event and may represent the final destiny of the X-ray binary SS 433.