Many interesting empirical models involve multiple discrete, often binary, choices or steps. Common examples include strategic, instrumental variable, and selection probits. Despite an abundance of binary time-series-cross-sectional data in political science, few applications with multiple equations account for unit-specific unobserved heterogeneity. One likely reason for this omission is that many of these models are numerically fickle and including large numbers of unit dummies can induce severe numerical instability. Indeed, multiple equation models compound standard incidental parameter problems, because each new equation comes with a full set of unit dummies. To overcome this problem, I propose a Mundlak specification to model unobserved heterogeneity in three advanced, but commonly used, models that are all based on systems of two binary equations. In all cases, the Mundlak specification outperforms dummy variable approaches even when the number of within-unit observations is quite large. The Mundlak approaches provide a simple and general solution for researchers to accommodate unit-specific heterogeneity without introducing increased numeric instability or excessive numbers of parameters, and they can be implemented using existing, commonly used software.