We have proposed new algorithms for the numerical integration of the equations of motion for classical spin systems. In close analogy to symplectic integrators for Hamiltonian equations of motion used in Molecular Dynamics these algorithms are based on the Suzuki-Trotter decomposition of exponential operators and unlike more commonly used algorithms exactly conserve spin length and, in special cases, energy. Using higher order decompositions we investigate integration schemes of up to fourth order and compare them to a well established fourth order predictor-corrector method. We demonstrate that these methods can be used with much larger time steps than the predictor-corrector method and thus may lead to a substantial speedup of computer simulations of the dynamical behavior of magnetic materials.