We present a method to treat the two-dimensional (2D) Hubbard model for parameter regimes which are relevant for the physics of the high-${T}_{c}$ superconducting cuprates. Unlike previous attempts to attack this problem, our approach takes into account all fluctuations in different channels on equal footing and is able to treat reasonable large lattice sizes up to $32\ifmmode\times\else\texttimes\fi{}32$. This is achieved by the following three-step procedure: (i) We transform the original problem to a representation (dual fermions) in which all purely local correlation effects from the dynamical mean field theory are already considered in the bare propagator and bare interaction of the new problem. (ii) The strong $1/{(i\ensuremath{\nu})}^{2}$ decay of the bare propagator allows us to integrate out all higher Matsubara frequencies aside from the lowest using low-order diagrams. The new effective action depends only on the two lowest Matsubara frequencies which allows us to (iii) apply the two-particle self-consistent parquet formalism, which takes into account the competition between different low-energy bosonic modes in an unbiased way, on much finer momentum grids than usual. In this way, we were able to map out the phase diagram of the 2D Hubbard model as a function of temperature and doping. Consistently with the experimental evidence for hole-doped cuprates and previous dynamical cluster approximation calculations, we find an antiferromagnetic region at low doping and a superconducting dome at higher doping. Our results also support the role of the van Hove singularity as an important ingredient for the high value of ${T}_{c}$ at optimal doping. At small doping, the destruction of antiferromagnetism is accompanied by an increase of charge fluctuations supporting the scenario of a phase-separated state driven by quantum critical fluctuations.