  We study the stability of neutron stars with toroidal magnetic fields by
magnetohydrodynamic simulation in full general relativity under assumption of
axial symmetry. Nonrotating and rigidly rotating neutron stars are prepared for
a variety of magnetic field configuration. For modeling the neutron stars, the
polytropic equation of state with the adiabatic index $\Gamma=2$ is used for
simplicity. It is found that nonrotating neutron stars are dynamically unstable
for the case that toroidal magnetic field strength varies $\propto
\varpi^{2k-1}$ with $k\geq 2$ (here $\varpi$ is the cylindrical radius),
whereas for $k=1$ the neutron stars are stable. After the onset of the
instability, unstable modes grow approximately in the Alfv\'en time scale and,
as a result, a convective motion is excited to change the magnetic field
profile until a new state, which is stable against axisymmetric perturbation,
is reached. We also find that rotation plays a role in stabilization, although
the instability still sets in in the Alfv\'en time scale when the ratio of
magnetic energy to rotational kinetic energy is larger than a critical value
$\sim 0.2$. Implication for the evolution of magnetized protoneutron stars is
discussed.
