  Self-interacting dark matter (SIDM) theory predicts that dark matter halos
experience core-collapse, a process where the halo's inner region rapidly
increases in density and decreases in size. The N-body simulations used to
study this process can suffer from numerical errors when simulation parameters
are selected incorrectly. Optimal choices for simulation parameters are well
studied for cold dark matter (CDM), but are not deeply understood when
self-interactions are included. In order to perform reliable N-body simulations
and model core-collapse accurately we must understand the potential numerical
errors, how to diagnose them, and what parameter selections must be made to
reduce them. We use the \texttt{Arepo} N-body code to perform convergence tests
of core-collapsing SIDM halos across a range of halo concentrations and SIDM
cross-sections, and quantify potential numerical issues related to mass
resolution, timestep size, and gravitational softening length. Our tests
discover that halos with fewer than $10^5$ simulation particles, a resolution
typically not met by subhalos in N-body simulations, suffer from significant
discreteness noise that leads to variation and extreme outliers in the collapse
rate. At our lowest resolution of $N=10^4$ particles, this collapse time
variation can reach as high as 20\%. At this low resolution we also find a bias
in collapse times and a small number of extreme outliers. Additionally, we find
that simulations which run far beyond the age of the Universe, which have been
used to calibrate SIDM gravothermal fluid models in previous work, have a
sensitivity to the timestep size that is not present in shorter simulations or
simulations using only CDM. Our work shows that choices of simulation
parameters that yield converged results for some halo masses and SIDM models do
not necessarily yield convergence for others.
