  We provide an algorithm for evolving general spin-$s$ Gross-Pitaevskii /
non-linear Schr\"odinger systems carrying a variety of interactions, where the
$2s+1$ components of the `spinor' field represent the different
spin-multiplicity states. We consider many nonrelativistic interactions up to
quartic order in the Schr\"odinger field (both short and long-range, and
spin-dependent and spin-independent interactions), including explicit
spin-orbit couplings. The algorithm allows for spatially varying external
and/or self-generated vector potentials that couple to the spin density of the
field. Our work can be used for scenarios ranging from laboratory systems such
as spinor Bose-Einstein condensates (BECs), to cosmological/astrophysical
systems such as self-interacting bosonic dark matter. As examples, we provide
results for two different setups of spin-$1$ BECs that employ a varying
magnetic field and spin-orbit coupling, respectively, and also collisions of
spin-$1$ solitons in dark matter. Our symplectic algorithm is second-order
accurate in time, and is extensible to the known higher-order accurate methods.
