  A systematically improvable variational wave function is proposed for the
numerical solution of strongly correlated systems, such as the celebrated
Hubbard model. With a stochastic optimization method, based on recent advances
in variational Monte Carlo, it is possible to achieve unbiased results in the
thermodynamic limit with simple and stable low order polynomial extrapolations
to zero temperature.This is achieved much before the so called sign problem
gets prohibitive. The convergence to the thermodynamic limit is indeed shown to
be possible already for relatively small square lattices with linear dimension
$L \simeq 16$, thanks to appropriate averages over several twisted boundary
conditions. Results for the phase diagram of the Hubbard model at $U/t=8$ are
reported. As a function of the doping from the half-filled Mott insulating
state, several fractional dopings $\delta =1/8,1/7,1/6$ are identified where
stripes of length $\xi ={1\over \delta}$ have insulating character. Then phase
separation is found to connect these antiferromagnetic phases toghether with
the half filled commensurate one. Moreover, the d wave uniform state, though
very close to the ground state energy, is clearly stable only in the large
doping region $>19\%$.
