Using numerical methods we investigate the dynamical stability of the Gliese 581 exoplanetary system. The system is known to harbour four planets (b-e). The existence of another planet (g) in the liquid water habitable zone of the star is debated after the latest analyses of the radial velocity (RV) measurements. We integrated the 4 and 5-planet model of Vogt et al. (AN 333, 561-575, 2012) with initial circular orbits. To characterize stability, the maximum eccentricity was used that the planets reached over the time of the integrations and the LCI and RLI to identify chaotic motion. Since circular orbits in the RV fits seem to be a too strong restriction and the true orbits might be elliptic, we investigated the stability of the planets as a function of their eccentricity. The integration of the circular 4-planet model shows that it is stable on a longer timescale for even an inclination i = 5{\deg}. A fifth planetary body in the 4-planet model could have a stable orbit between the two super-Earth sized planets c and d, and beyond the orbit of planet d, although another planet would likely only be stable on circular or near-circular orbit in the habitable zone of the star. Gliese 581 g in the 5-planet model would have a dynamically stable orbit, even for a wider range of orbital parameters, but its stability is strongly dependent on the eccentricity of planet d. The low-mass planet e, which quickly became unstable in eccentric models, remains stable in the circular 4-planet model, but the stable region around its initial semi-major axis and eccentricity is rather small. The stability of the inner planets e and c is dependent on the eccentricity of the Neptune-size planet b. The outermost planet d is far away from the adjacent planet c to considerably influence its stability, however, the existence of a planet between the two super-Earth planets c and d constrains its eccentricity.