Yes mean-field theory is wrong for the one-dimensional case (and wrong for the two and three dimensional cases as well, where the transition exists but the mean-field approximation gets the wrong critical temperature and exponents). In fact it's a typical first year exercise to solve the 1D Ising model exactly using transfer matrices, and I suggest you look into that.
The nature of the mean-field approximation is that it assumes there are no thermal fluctuations around the approximate solution you propose (i.e., a state with ferromagnetic order) but in low dimensions, this approximation is often qualitatively wrong.
The mean-field theory of the Ising model happens to be exact in 4-dimensions, but more complicated phase transitions might not be well described by mean-field theory for even higher dimensions (this is called the "upper critical dimension").
I will explain how I measured the spin-spin correlation function for the 2d Ising model. Generalization to more than 2 dimensions should be straightforward as long as you have hypercubic lattices.
Just to get the notation straight: Let's use the name $\sigma_{(i,j)}$ for the spin at position $\vec r_{(i,j)} = \begin{pmatrix} i \\ j \end{pmatrix}$.
Let's assume our Ising model has $L \times L$ spins, so we have $1 \leq i,j \leq L$.
The first thing to note is that I assume that the spin-spin correlation function - lets call it $\chi_{(i,j)(k,l)}$ - only depends on the absolute distance between the two spins $\sigma_{(i,j)}$ and $\sigma_{(k,l)}$.
$$\chi_{(i,j)(k,l)} = \chi( | \vec r_{(i,j)} - \vec r_{(k,l)} | ) = \chi(r)$$
I'm not sure if this is true in general, but there is a proof that it is true for the 2d Ising model at the critical temperature. (A rotationally symmetric spin-spin correlation function is proven, and that is just another way of saying that it depends on the distance only.) Don't take my word for it, but for now I wouldn't worry about this assumption ...
Knowing this, it is easy to see that at every step in your Markov chain you don't have to look at all possible $L^2$ pairs of spins. Instead you can only look at a subset of pairs, where you have all distances $r$ you are interested in within this subset.
For my implementation I chose to use all spins $\sigma_{(i,i)}$ along the diagonal of the lattice and to only measure the correlation along the two directions of the lattice vectors. So I only consider the correlation between $\sigma_{(i,i)}$ and all $\sigma_{(i,j)}$ (horizontally) and $\sigma_{(j,i)}$ (vertically). Remember that the spin-spin correlation function is rotationally symmetric, so you don't gain anything by considering pairs of spins along any other direction. Measuring along the axes makes your life much easier as you will only have integer distances $r$, meaning that you can use a plain $L-1$ element array with the distance as an index to accumulate the terms of your sum. This is what you should calculate at every Monte Carlo step:
$$ \text{sum}(r) = \sum_{i=1}^L \sum_{j=1}^L \left( \sigma_{(i,i)} \sigma_{(i,j)} + \sigma_{(i,i)} \sigma_{(j,i)} \right) \; \text{ with } \; r = |i-j|$$
$$ \text{samples}(r) = \sum_{i=1}^L \sum_{j=1}^L 2 \; \text{ with } \; r = |i-j|$$
$$ \chi_n(r) = \frac{\text{sum}(r)}{\text{samples}(r)}$$
This is exactly what this function IsingModel2d::ss_corr() in my code does. (Keeping track of the samples is only needed because you make a different number of measurements for different values of $r$. Consider as an example the extreme cases: All but two spins $\sigma_{(i,i)}$ on the diagonal have 4 nearest neighbors at distance $r=1$, but only $\sigma_{(1,1)}$ and $\sigma_{(L,L)}$ have two partners each at distance $r = L-1$. So you have to compensate for the fact that you measure the correlation between close spins more often than between far away spins.)
Note that you will get a different $\chi_n(r)$ at every Monte Carlo step, depending on your current configuration. In the end you should average them (just like you do with any other observable) to get your final result for the spin-spin correlation function:
$$\chi(r) = \frac{1}{N-1} \sum_{n=1}^N \chi_n(r)$$
Once you have this, you can fit $\chi(r) = r^{-\beta}$ where you should hopefully find $\beta = 0.25$. Check out this other question on details regarding the fit.
Hope this helps!
Best Answer
That is not a definition of correlation length. (It is a definition of the critical exponent.)
The correlation length is defined in terms of the 2-point correlation function of spin observables. Pick points $x$ and $y$ on the lattice, and consider the expectation value $\langle s(x) s(y) \rangle$ of the product of the spin observable at $x$ and the spin observable at $y$. This quantity tells you how strongly correlated the spin at $x$ and the spin at $y$ are, as a function of the temperature, coupling constant, and the distance between $x$ and $y$. If $T > T_c$, then the correlation function dies off exponentially fast in $|x-y|$.
$\langle s(x) s(y) \rangle \sim e^{-\frac{|x-y|}{\xi(T)}}$
The correlation length is, by definition, the constant (in $x$ and $y$, but not in $T$) which tells you how fast the correlation function vanishes.