Variation in a sample of molecular sequence data informs about the past evolutionary history of the sample's population. Traditionally, Bayesian modelling coupled with the standard coalescent is used to infer the sample's bifurcating genealogy and demographic and evolutionary parameters such as effective population size and mutation rates. However, there are many situations where binary coalescent models do not accurately reflect the true underlying ancestral processes. Here, we propose a Bayesian non-parametric method for inferring effective population size trajectories from a multifurcating genealogy under the [Formula: see text]-coalescent. In particular, we jointly estimate the effective population size and the model parameter for the Beta-coalescent model, a special type of [Formula: see text]-coalescent. Finally, we test our methods on simulations and apply them to study various viral dynamics as well as Japanese sardine population size changes over time. The code and vignettes can be found in the phylodyn package.This article is part of the theme issue '"A mathematical theory of evolution": phylogenetic models dating back 100 years'.