We revisit Hawkins’ (Comput Stat 9(3):233–247, <CitationRef CitationID="CR15">1994</CitationRef>) algorithm for fitting monotonic polynomials and discuss some practical issues that we encountered using this algorithm, for example when fitting high degree polynomials or situations with a sparse design matrix but multiple observations per <InlineEquation ID="IEq1"> <EquationSource Format="TEX">$$x$$</EquationSource> </InlineEquation>-value. As an alternative, we describe a new approach to fitting monotone polynomials to data, based on different characterisations of monotone polynomials and using a Levenberg–Marquardt type algorithm. We consider different parameterisations, examine effective starting values for the non-linear algorithms, and discuss some limitations. We illustrate our methodology with examples of simulated and real world data. All algorithms discussed in this paper are available in the R Development Core Team (A language and environment for statistical computing, R Foundation for Statistical Computing, Vienna, <CitationRef CitationID="CR25">2011</CitationRef>) package <Emphasis FontCategory="NonProportional">MonoPoly. Copyright Springer-Verlag Berlin Heidelberg 2013
