The long-studied Hubbard model is one of the simplest models of copper-oxide superconductors. However, the connection between the model and the experimental phase diagram is still under debate, in particular regarding the existence and extent of the $d$-wave superconducting phase. Recent rapid progress in improving the accuracy of numerical solvers has opened a way to answer this question reliably. Here, we study the hole-doping concentration ($\delta$) dependence of the Hubbard model in the ground states on a square lattice at strong coupling $U/t=10$, for the on-site interaction $U$ and the transfer $t$, using a variational Monte Carlo method. The method, which combines tensor network and Lanczos methods on top of Pfaffian wave functions, reveals a rich phase diagram, in which many orders compete severely and degenerate within the energy range of 0.01$t$. We have identified distinct phases including a uniform $d$-wave superconducting phase for $0.17\lesssim \delta \lesssim0.22$ and a stripe charge/spin ordered phase for $\delta\lesssim0.17$ with the stripe period depending on $\delta$, together with presumable spatially coexisting antiferromagnetic and stripe order for $\delta\lesssim0.07$ and coexisting stripe and $d$-wave superconductivity for $0.07\lesssim\delta\lesssim0.17$. The present, improved method revealed a wider region of a charge uniform superconducting phase than the previous studies and shows a qualitative similarity to the phase diagram of the cuprate superconductors. The superconducting order parameter is largest at doping of around $\delta=0.17$ in the ground state, which undergoes phase transitions from an inhomogeneous to a uniform state.