<p style='text-indent:20px;'>The dynamics of many systems from physics, economics, chemistry, and biology can be modelled through polynomial functions. In this paper, we provide a computational means to find positively invariant sets of polynomial dynamical systems by using semidefinite programming to solve <i>sum-of-squares (SOS)</i> programmes. With the emergence of SOS programmes, it is possible to efficiently search for Lyapunov functions that guarantee stability of polynomial systems. Yet, SOS computations often fail to find functions, such that the conditions hold in the entire state space. We show here that restricting the SOS optimisation to specific domains enables us to obtain positively invariant sets, thus facilitating the analysis of the dynamics by considering separately each positively invariant set. In addition, we go beyond classical Lyapunov stability analysis and use SOS decompositions to computationally implement sufficient positivity conditions that guarantee existence, uniqueness, and exponential stability of a limit cycle. Importantly, this approach is applicable to systems of any dimension and, thus, goes beyond classical methods that are restricted to two-dimensional phase space. We illustrate our different results with applications to classical systems, such as the van der Pol oscillator, the Fitzhugh-Nagumo neuronal equation, and the Lorenz system.</p>