The reachable set of controlled dynamical systems is the set of all reachable states from an initial condition over a certain time horizon, subject to operational constraints and exogenous disturbances. In astrodynamics, rapid approximation of reachable sets is invaluable for trajectory planning, collision avoidance, and ensuring safe and optimal performance in complex dynamics. Leveraging the connection between minimum-time trajectories and the boundary of reachable sets, we propose a sampling-based method for rapid and efficient approximation of reachable sets for finite- and low-thrust spacecraft. The proposed method combines a minimum-time multi-stage indirect formulation with the celebrated primer vector theory. Reachable sets are generated under two-body and circular restricted three-body (CR3B) dynamics. For the two-body dynamics, reachable sets are generated for (1) the heliocentric phase of a benchmark Earth-to-Mars problem, (2) two scenarios with uncertainties in the initial position and velocity of the spacecraft at the time of departure from Earth, and (3) a scenario with a bounded single impulse at the time of departure from Earth. For the CR3B dynamics, several cislunar applications are considered, including L1 Halo orbit, L2 Halo orbit, and Lunar Gateway 9:2 NRHO. The results indicate that low-thrust spacecraft reachable sets coincide with invariant manifolds existing in multi-body dynamical environments. The proposed method serves as a valuable tool for qualitatively analyzing the evolution of reachable sets under complex dynamics, which would otherwise be either incoherent with existing grid-based reachability approaches or computationally intractable with a complete Hamilton–Jacobi–Bellman method.