We employ the nonperturbative many-body spin GW method to compute the single-particle Green's function, and the covariance method to evaluate the two-body correlation functions, for the Hubbard model of superconductors. The method, benchmarked at strong coupling and low temperatures, allows simulations on lattices up to
64×64, thereby enabling access to the thermodynamic limit. We find the key transition temperatures---namely the Néel temperature and the superconducting critical temperature
Tc---to be in quantitative agreement with experimental data. Above the
d-wave superconducting dome, there exists a broad pseudogap region in which the normal-state excitations are suppressed. These results indicate our approach may offer a new route toward investigating strongly correlated systems, especially cuprate superconductors.