Spatially-resolved stellar populations are the gold-standard for studying star formation and stellar evolution, from star clusters in the Milky Way to galaxies in the Local Group. By measuring the photometry of individual stars and constructing a colour-magnitude diagram (CMD), we gain exquisite resolution into the stellar masses, ages, and metallicities present in a population. This enables us to probe fundamental questions about the stellar birth function (SBF), a complete historical account of the star formation in an observed population. The upcoming James Webb Space Telescope (JWST) and Roman Space Telescope (Roman) will provide an unprecedented photometric depth and wide field of view respectively, increasing the number of observed stars from any given population by several orders of magnitude. There is an urgent need for new data analysis techniques that can take advantage of these forthcoming massive datasets. We present our ongoing development of a suite of software tools to fit the stellar birth function (including the initial mass function, star formation history, and metallicity distribution) of resolved stellar populations. We model the underlying stellar population as a Poisson point process from which stars on the CMD are drawn, and compare the observed population to forward-modelled theoretical populations with a novel distance metric. We circumvent the pitfalls of traditional binning and exact-counting methods by employing an approximate Bayesian algorithm to derive the posterior distributions of stellar birth function parameters without ever defining an explicit likelihood. Early simulation studies demonstrate that our method improves upon the accuracy of past techniques whilst using a fraction of the computational resources, especially in the limit of many stars on the CMD. We are currently using our software to analyze the Small Magellanic Cloud using Hubble Space Telescope data. We will release and distribute our package in time for JWST early release science programs, allowing the community to test it on globular clusters (M92), ultra-faint dwarfs (Draco II), and star-forming dwarfs (WLM) with JWST data.