We present the first simulated galaxy clusters (M200 > 1014 M⊙) with both self-interacting dark matter (SIDM) and baryonic physics. They exhibit a greater diversity in both dark matter and stellar density profiles than their counterparts in simulations with collisionless dark matter (CDM), which is generated by the complex interplay between dark matter self-interactions and baryonic physics. Despite variations in formation history, we demonstrate that analytical Jeans modelling predicts the SIDM density profiles remarkably well, and the diverse properties of the haloes can be understood in terms of their different final baryon distributions.