We calculate axisymmetric toroidal modes of magnetized neutron stars with a solid crust. We assume the interior of the star is threaded by a poloidal magnetic field that is continuous at the surface with the outside dipole field whose strength Bp at the magnetic pole is Bp ∼ 10 16G. Since separation of variables is not possible for oscillations of magnetized stars, to represent the angular dependence of the oscillation modes we employ finite series expansions in terms of spherical harmonic functions. For Bp ∼ 1016G, we find isolated sequences of low radial order toroidal modes, whose oscillation frequency slowly increases as the number of radial nodes of the eigenfunction increases. The frequency spectrum of toroidal modes for Bp ∼ 1016G is largely different from that of crustal toroidal modes of the non-magnetized model.