  Fixed-node Green's function Monte Carlo calculations have been performed for
very large 16x6 2D Hubbard lattices, large interaction strengths U=10,20, and
40, and many (15-20) densities between empty and half filling. The nodes were
fixed by a simple Slater-Gutzwiller trial wavefunction. For each value of U we
obtained a sequence of ground-state energies which is consistent with the
possibility of a phase separation close to half-filling, with a hole density in
the hole-rich phase which is a decreasing function of U. The energies suffer,
however, from a fixed-node bias: more accurate nodes are needed to confirm this
picture. Our extensive numerical results and their test against size, shell,
shape and boundary condition effects also suggest that phase separation is
quite a delicate issue, on which simulations based on smaller lattices than
considered here are unlikely to give reliable predictions.
