We study the stability of neutron stars with toroidal magnetic fields by magnetohydrodynamic simulation in full general relativity under the 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 where toroidal magnetic field strength varies proportional to pi(2k-1) with k >= 2 (here pi 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 Alfven 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 occurs in the Alfven time scale when the ratio of magnetic energy to rotational kinetic energy is larger than a critical value similar to 0.2. Implication for the evolution of magnetized protoneutron stars is discussed.