AbstractContext: The M dwarf Gliese 581 has recently been found to harbour two super Earths in addition to an already known close-in Neptune-mass planet. Interestingly, these two planets are considered as potentially habitable, and recent theoretical works give further credit to this hypothesis, in particular for the outermost planet(Gl 581 d). Aims: In this paper, we address the issue of the dynamical stability evolution of this planetary system. This is important because the basic stability ensures that a 3-planet model is a physically adequate description of the radial-velocity(RV) data. It is also crucial when considering the planets' habitability because the secular evolution of the orbits may regulate their climate, even in the case where the system is stable. Methods: We have numerically integrated the planetary system over 108 yr, starting from the present fitted solution. We also performed additional simulations where i) we vary the inclination of the system relative to the line of sight, ii) assume eccentricities at the upper limit of the error bars in the radial velocity fit and where iii) we consider additional(yet undetected) outer planets. We also compute Lyapunov exponents to quantify the level of dynamical chaos in the system. Results: In all cases, the system appears dynamically stable, even in close to pole-on configurations. The system is actually chaotic, but nevertheless stable. The semi-major axes of the planets are extremely stable, and their eccentricities undergo small amplitude variations. The addition of potential outer planets does not affect this result. Conclusions: Consequently, from the dynamical point-of-view, a 3-planet model is an adequate description of the present RV-data set. Only a limited range of inclinations can be excluded for coplanar orbits(i<10░). The climate on the planets is expected to be secularly stable, thus not precluding the development of life. Gl 581 remains the best candidate for a planetary system with planets that potentially bear primitive forms of life.