  We study numerically the two-point correlation functions of height functions
in the six-vertex model with domain wall boundary conditions. The correlation
functions and the height functions are computed by the Markov chain Monte-Carlo
algorithm. Particular attention is paid to the free fermionic point
($\Delta=0$), for which the correlation functions are obtained analytically in
the thermodynamic limit. A good agreement of the exact and numerical results
for the free fermionic point allows us to extend calculations to the disordered
($|\Delta|<1$) phase and to monitor the logarithm-like behavior of correlation
functions there. For the antiferroelectric ($\Delta<-1$) phase, the exponential
decrease of correlation functions is observed.
