  When dark matter has a large cross section for self scattering, halos can
undergo a process known as gravothermal core collapse, where the inner core
rapidly increases in density and temperature. To date, several methods have
been used to implement Self-Interacting Dark Matter~(SIDM) in N-body codes, but
there has been no systematic study of these different methods or their accuracy
in the core-collapse phase. In this paper, we compare three different numerical
implementations of SIDM, including the standard methods from the GIZMO and
Arepo codes, by simulating idealized dwarf halos undergoing significant dark
matter self interactions ($\sigma/m = 50$~cm$^2$/g). When simulating these
halos, we also vary the mass resolution, time-stepping criteria, and
gravitational force-softening scheme. The various SIDM methods lead to distinct
differences in a halo's evolution during the core-collapse phase, as each
results in slightly different scattering rates and spurious energy
gains/losses. The use of adaptive force softening for gravity can lead to
numerical heating that artificially accelerates core collapse, while an
insufficiently small simulation time step can cause core evolution to stall or
completely reverse. Additionally, particle numbers must be large enough to
ensure that the simulated halos are not sensitive to noise in the initial
conditions. Even for the highest-resolution simulations tested in this study
($10^6$ particles per halo), we find that variations of order $10\%$ in
collapse time are still present. The results of this work underscore the
sensitivity of SIDM modeling on the choice of numerical implementation and
motivate a careful study of how these results generalize to halos in a
cosmological context.
