In the definition of a positive-definite function, why are the elements of the n$\times$n matrix chosen as $f(x_i-x_j)$ for $i, j = 1,...,n$?
Also, it says that "for any real numbers $x_i$". Does this mean that this works only for functions that are defined over the entire $\mathbb{R}$?
The reason for considering $f(x_i-x_j)$ is stated just below on the same page: Bochner's theorem, which characterizes such functions as Fourier transforms of positive measures. This is the result that motivates the definition.
The kernel of the Fourier transform, namely the exponential function $e(x,t)=e^{-2\pi ix t}$, has the property that $$ e(x_i-x_j,t) = e(x_i ,t)\overline{e(x_j,t)} $$ Hence, the matrix with entries $e(x_i-x_j,t)$ is of the form $vv^*$ where $ v $ is the column vector with entries $e(x_i,t)$. This implies that the matrix is positive semidefinite.
Since positive semidefinite matrices form a convex cone, one half of Bochner's theorem readily follows: the Fourier transform of any positive measure is a positive-definite function. (The converse is the hard part of the theorem.)
The concept is not a natural one for functions that are defined only on a subinterval of $\mathbb{R}$. Its root is in the additive group structure of $\mathbb{R}$, which finite intervals do not have. It makes sense to consider it for functions on additive subgroups of $\mathbb{R}$, such as $\mathbb{Z}$. Or even on more general groups.