We present methods to fit monotone polynomials in a Bayesian framework, including implementations in the popular, readily available, modeling languages BUGS and Stan. The sum-of-squared polynomials parameterisation of monotone polynomials previously considered in the frequentist framework by Murray, Müller & Turlach (2016), is again considered here, due to its superior flexibility compared to other parameterisations. The specifics of our implementation are discussed, enabling end users to adapt this work to their applications. Testing was undertaken on real and simulated data sets, the output and diagnostics of which are presented. We demonstrate that Stan is preferable for high degree polynomials, with the component-wise nature of Gibbs sampling being potentially inappropriate for such highly connected models. All code discussed here, and sample scripts that show how to use it from R, is freely available at https://github.com/hhau/BayesianMonPol.