In this work, we highlight how trapped-ion quantum systems can be used to study generalized Holstein models, and benchmark expensive numerical calculations. We study a particular spin-Holstein model that can be implemented with arrays of ions confined by individual microtraps, and that is closely related to the Holstein model of condensed matter physics, used to describe electron-phonon interactions. In contrast to earlier proposals, we focus on realizing many-electron systems and inspect the competition between charge-density wave order, fermion pairing and phase separation. In our numerical study, we employ a combination of complementary approaches, based on non-Gaussian variational ansatz states and matrix product states, respectively. We demonstrate that this hybrid approach outperforms standard density-matrix renormalization group calculations.