There is a nice review Numerical studies of the 2D Hubbard model by Doug Scalapino which considers the question.
He claims that in the end the answer is yes.
However, there are several subtleties along the way. I just mention two that I learnt.
To attempt to answer the question one calculates the d-wave pairing susceptibility on a finite lattice.
the pairing susceptibility is less than for U=0!
Clearly there is no superconductivity for U=0 and a key claim of the Anderson paradigm is that large U and correlations are the key for high-Tc. So what is going on?
Scalapino and his collaborators argue that this decrease with increasing U is because U reduces the single particle spectral weight [i.e. the effective mass of quasi-particles increases with increasing U] and so P_d should be compared to the function below which is the value of P_d in the absence of interactions between the holes.
finite size effects and dependence on cluster geometry
It turns out that the results don't just depend on the size of the cluster but also its geometry. "One needs to take into account the 4-site local plaquette structure of the order parameter." Z_d is the number of independent near-neighbour plaquettes on a cluster. For the square 16 site cluster (16B) in the figure below Z_d=2. In contrast, for the diamond 16 site cluster (16A) on has Z_d=3. As a result the pairing correlations are suppressed in 16A relative to 16B.
Getting or not getting superconductivity in quantum Monte simulations is certainly subtle. A different technique [constrained path Monte Carlo and which is at zero temperature] used by a different group did not produce long range correlations. [See this PRB from 1999]. I am not sure if this discrepancy has ever been explained to the satisfaction of all parties...