We probe the self-interactions of dark matter using observational data of relaxed galaxy groups and clusters. Our analysis uses the Jeans formalism and considers a wider range of systematic effects than in previous work, including adiabatic contraction and stellar anisotropy, to robustly constrain the self-interaction cross section. For both groups and clusters, our results show a mild preference for a nonzero cross section compared with cold collisionless dark matter. Our groups result, $\sigma/m=0.5\pm0.2~\mathrm{cm}^2/\mathrm{g}$, places the first constraint on self-interacting dark matter (SIDM) at an intermediate scale between galaxies and massive clusters. Our clusters result is $\sigma/m=0.19\pm0.09~\mathrm{cm}^2/\mathrm{g}$, with an upper limit of $\sigma / m < 0.35~\mathrm{cm}^2/\mathrm{g}$ (95% CL). Thus, our results disfavor a velocity-independent cross section of order $1~\mathrm{cm}^2/\mathrm{g}$ or larger needed to address small scale structure problems in galaxies, but are consistent with a velocity-dependent cross section that decreases with increasing scattering velocity. Comparing the cross sections with and without the effect of adiabatic contraction, we find that adiabatic contraction produces slightly larger values for our data sample, but they are consistent at the $1\sigma$ level. Finally, to validate our approach, we apply our Jeans analysis to a sample of mock data generated from SIDM-plus-baryons simulations with $\sigma/m = 1~\mathrm{cm}^2/\mathrm{g}$. This is the first test of the Jeans model at the level of stellar and lensing observables directly measured from simulations. We find our analysis gives a robust determination of the cross section, as well as consistently inferring the true baryon and dark matter density profiles.