We derive second and higher order explicit symplectic integrators for the charged particle motion in an s-dependent magnetic field with the paraxial approximation. The Hamiltonian of such a system takes the form of $H = \sum_k(p_k-a_k(\vec{q},s))^2 + V(\vec{q},s)$. This work solves a long-standing problem for modeling $s$-dependent magnetic elements. Important applications of this work include the studies of the charged particle dynamics in a storage ring with strong field wigglers, arbitrarily polarized insertion devices,and super-conducting magnets with strong fringe fields. Consequently, this work will have a significant impact on the optimal use of the above magnetic devices in the light source rings as well as in next generation linear collider damping rings.