  Maximally Smooth Functions (MSFs) are a form of constrained functions in
which there are no inflection points or zero crossings in high order
derivatives. Consequently, they have applications to signal recovery in
experiments where signals of interest are expected to be non-smooth features
masked by larger smooth signals or foregrounds. They can also act as a powerful
tool for diagnosing the presence of systematics. The constrained nature of MSFs
makes fitting these functions a non-trivial task. We introduce maxsmooth, an
open source package that uses quadratic programming to rapidly fit MSFs. We
demonstrate the efficiency and reliability of maxsmooth by comparison to
commonly used fitting routines and show that we can reduce the fitting time by
approximately two orders of magnitude. We introduce and implement with
maxsmooth Partially Smooth Functions, which are useful for describing elements
of non-smooth structure in foregrounds. This work has been motivated by the
problem of foreground modelling in 21-cm cosmology. We discuss applications of
maxsmooth to 21-cm cosmology and highlight this with examples using data from
the Experiment to Detect the Global Epoch of Reionization Signature (EDGES) and
the Large-aperture Experiment to Detect the Dark Ages (LEDA) experiments. We
demonstrate the presence of a sinusoidal systematic in the EDGES data with a
log-evidence difference of $86.19\pm0.12$ when compared to a pure foreground
fit. MSFs are applied to data from LEDA for the first time in this paper and we
identify the presence of sinusoidal systematics. maxsmooth is pip installable
and available for download at: https://github.com/htjb/maxsmooth
