Bell non-locality represents the ultimate consequence of quantum entanglement, fundamentally undermining the classical tenet that spatially-separated degrees of freedom possess objective attributes independently of the act of their measurement. Despite its importance, probing Bell non-locality in many-body systems is considered to be a formidable challenge, with a computational cost scaling exponentially with system size. Here we propose and validate an efficient variational scheme, based on the solution of inverse classical Ising problems, which in polynomial time can probe whether an arbitrary set of quantum data is compatible with a local theory; and, if not, it delivers a many-body Bell inequality violated by the quantum data. We use our approach to unveil new many-body Bell inequalities, violated by suitable measurements on paradigmatic quantum states (the low-energy states of Heisenberg antiferromagnets), paving the way to systematic Bell tests in the many-body realm.