Abstract: in order to prepare for the huge amount of data incoming from future 21cm missions, we develop a fully analytical approach to studying the ionized bubble size distribution (BSD) during the epoch of reionization, motivated by seminal work by Furlanetto at al (2004). We develop this as part of a package known as Zeus21 that interfaces with CLASS cosmology to compute the 21cm spin temperature and global signal. The addition of reionization bubbles completes the calculation by adding fluctuations into the signal. Using our technique, we provide a way to scan large portions of parameter space in a few seconds, complementing detailed large-scale hydrodynamic simulations.