Unveiling the evolution of toroidal field instability, known as Tayler instability, is essential to understand the strength and topology of the magnetic fields observed in early-type stars, in the core of the red giants, or in any stellar radiative zone. We want to study the non-linear evolution of the instability of a toroidal field stored in a stably stratified layer, in spherical symmetry and in the absence of rotation. In particular, we intend to quantify the suppression of the instability as a function of the Brunt-Väisäla (ω BV ) and the Alfvén (ω A ) frequencies. We use the MHD equations as implemented in the anelastic approximation in the EULAG-MHD code and perform a large series of numerical simulations of the instability exploring the parameter space for the ω BV and ω A . We show that beyond a critical value gravity strongly suppress the instability, in agreement with the linear analysis. The intensity of the initial field also plays an important role: weaker fields show much slower growth rates. Moreover, in the case of very low gravity, the fastest growing modes have a large characteristic radial scale, at variance with the case of strong gravity, where the instability is characterized by horizontal displacements. Our results illustrate that the anelastic approximation can efficiently describe the evolution of toroidal field instability in stellar interiors. The suppression of the instability as a consequence of increasing values of ω BV might play a role to explain the magnetic desert in Ap/Bp stars, since weak fields are only marginally unstable in the case of strong gravity.