We introduce an approach to find approximate numerical solutions of truncated bootstrap equations for conformal field theories (CFTs) in arbitrary dimensions. The method is based on a stochastic search via a Metropolis algorithm guided by an action S which is the logarithm of the truncated bootstrap equations for a single scalar field correlator. While numerical conformal bootstrap methods based on semidefinite programming put rigorous exclusion bounds on CFTs, this method looks for approximate solutions, which correspond to local minima of S, when present, and can be even far from the extremality region. By this protocol we find that if no constraint on the operator scaling dimensions is imposed, S has a single minimum, corresponding to the free theory. If we fix the external operator dimension, however, we encounter minima that can be studied with our approach. Imposing a conserved stress-tensor, a Z(2) symmetry and one relevant scalar, we identify two regions where local minima of S are present. When projected in the (Delta(sigma), Delta(epsilon))-plane, sigma and epsilon being the external and the lightest exchanged operators, one of these regions essentially coincides with the extremality line found in previous bootstrap studies. The other region is along the generalized free theories in d = 2 and below that in both d = 3 and d = 4. We empirically prove that some of the minima found are associated to known theories, including the 2d and 3d Ising theories and the 2d Yang-Lee model.
