We present a new random approximation method that yields the existence of a discrete Beurling prime system P = {p 1 , p 2 , . . . } which is very close in a certain precise sense to a given non-decreasing, right-continuous, nonnegative, and unbounded function F . This discretization procedure improves an earlier discrete random approximation method due to H. Diamond, H. Montgomery, and U. Vorhauer [Math. Ann. 334 (2006), 1-36], and refined by W.-B. Zhang [Math. Ann. 337 (2007), 671-704].We obtain several applications. Our new method is applied to a question posed by M. Balazard concerning Dirichlet series with a unique zero in their half plane of convergence, to construct examples of very well-behaved generalized number systems that solve a recent open question raised by T. Hilberdink and A. Neamah in [Int. J. Number Theory 16 05 (2020), 1005-1011], and to improve the main result from [Adv. Math. 370 (2020), Article 107240], where a Beurling prime system with regular primes but extremely irregular integers was constructed. 2020 Mathematics Subject Classification. 11M41, 11N80. Key words and phrases. Discrete random approximation; Diamond-Montgomery-Vorhauer-Zhang probabilistic method; Dirichlet series with unique zero in half plane of convergence; Well-behaved Beurling primes and integers; Beurling integers with large oscillation; Riemann hypothesis for Beurling numbers.F.