Exospheres, the tenuous gas environments surrounding planets, planetary satellites, and cometary comae, play a significant role in mediating the interactions of these astronomical bodies with their surrounding space environments. This paper presents a comprehensive review of both analytical and numerical methods employed in modeling exospheres. The paper explores analytical models, including the Chamberlain and Haser models, which have significantly contributed to our understanding of exospheres of planets, planetary satellites, and cometary comae. Despite their simplicity, these models provide baselines for more complex simulations. Numerical methods, particularly the Direct Simulation Monte Carlo (DSMC) method, have proven to be highly effective in capturing the detailed dynamics of exospheres under non-equilibrium conditions. The DSMC method’s capacity to incorporate a wide range of physical processes, such as particle collisions, chemical reactions, and surface interactions, makes it an indispensable tool in planetary science. The Adaptive Mesh Particle Simulator (AMPS), which employs the DSMC method, has demonstrated its versatility and effectiveness in simulating gases in planetary and satellite exospheres and dusty gas cometary comae. It provides a detailed characterization of the physical processes that govern these environments. Additionally, the multi-fluid model BATSRUS has been effective in modeling neutral gases in cometary comae, as discussed in the paper. The paper presents methodologies of exosphere modeling and illustrates them with specific examples, including the modeling of the Enceladus plume, the sodium exosphere of the Moon, the coma of comet 67P/Churyumov-Gerasimenko, and the hot oxygen corona of Mars and Venus.