The onset of strange baryons in the core of neutron stars and the consequent softening of the equation of state have been questioned for a long time. Controversial theoretical predictions about the maximum mass of neutron stars and the recent astrophysical observations are the grounds for the so-called hyperon puzzle. We attempt to give our contribution to the discussion by studying the general problem of the hyperon-nucleon interaction by means of Auxiliary Field Diffusion Monte Carlo calculations. We employ a phenomenological approach showing that a three-body hyperon-nucleon force provides the strong repulsive contribution needed to correctly describe the systematics of medium-light Lambda-hypernuclei [1,2]. The same potential has been used to determine the equation of state and the mass-radius relation of an infinite systems of neutrons and Lambda particles. We find that the three-body hyperon-nucleon force has a dramatic effect on the equation of state and the predicted maximum mass [3]. Our results suggest that more constraints on the nature of hyperon-neutron forces are needed before drawing any conclusion on the role played by hyperons in neutron stars. [1] D. Lonardoni, S. Gandolfi, F. Pederiva, Phys. Rev. C 87, 041303(R) (2013) [2] D. Lonardoni, F. Pederiva, S. Gandolfi, Phys. Rev. C 89, 014314 (2014) [3] D. Lonardoni, A. Lovato, S. Gandolfi, F. Pederiva, arXiv:1407.4448 [nucl-th] (2014), submitted to Phys. Rev. Lett.