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 (Δ = 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 (|Δ| < 1) phase and to monitor the logarithm-like behavior of correlation functions there. For the antiferroelectric (Δ < −1) phase, the exponential decrease of correlation functions is observed.