  The onset of hyperons 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 and recent astrophysical observations of
neutron stars are the grounds for the so-called hyperon puzzle. We calculate
the equation of state and the neutron star mass-radius relation of an infinite
systems of neutrons and $\Lambda$ particles by using the auxiliary field
diffusion Monte Carlo algorithm. We find that the three-body hyperon-nucleon
interaction plays a fundamental role in the softening of the equation of state
and for the consequent reduction of the predicted maximum mass. We have
considered two different models of three-body force that successfully describe
the binding energy of medium mass hypernuclei. Our results indicate that they
give dramatically different results on the maximum mass of neutron stars, not
necessarily incompatible with the recent observation of very massive neutron
stars. We conclude that stronger constraints on the hyperon-neutron force are
necessary in order to properly assess the role of hyperons in neutron stars.
