  In recent years, observations have shown that multiple-star systems such as
hierarchical triple and quadruple-star systems are common, especially among
massive stars. They are potential sources of interesting astrophysical
phenomena such as compact object mergers, leading to supernovae, and
gravitational wave events. However, many uncertainties remain in their often
complex evolution. Here, we present the population synthesis code Multiple
Stellar Evolution (MSE), designed to rapidly model the stellar, binary, and
dynamical evolution of multiple-star systems. MSE includes a number of new
features not present in previous population synthesis codes: (1) an arbitrary
number of stars, as long as the initial system is hierarchical, (2) dynamic
switching between secular and direct N-body integration for efficient
computation of the gravitational dynamics, (3) treatment of mass transfer in
eccentric orbits, which occurs commonly in multiple-star systems, (4) a simple
treatment of tidal, common-envelope, and mass transfer evolution in which the
accretor is a binary instead of a single star, (5) taking into account planets
within the stellar system, and (6) including gravitational perturbations from
passing field stars. MSE, written primarily in the C++ language, will be made
publicly available and has few prerequisites; a convenient Python interface is
provided. We give a detailed description of MSE and illustrate how to use the
code in practice. We demonstrate its operation in a number of examples.
