We study numerically the density profile in the six-vertex model with domain wall boundary conditions. Using a Monte Carlo algorithm originally proposed by Allison and Reshetikhin we numerically evaluate the inhomogeneous density profiles in the disordered and antiferromagnetic regimes where frozen corners appear. At the free fermion point we present an exact finite-size formula for the density on the horizontal edges that relies on the imaginary time transfer matrix approach. In all cases where exact analytic forms for the density and the arctic curves are known the numerical method shows perfect agreement with them. This also suggests the possibility of its use for accurate quantitative purposes