At exactly 100 degrees C water is in a complicated state. There are bubbles of steam roiling beneath the surface, popping and condensing or evaporating away. It looks complicated. The boiling point is a critical point of the system. Critical points are worth studying and there is a simple model that encapsulates some of the physics that occur there – the Ising Model. The simulation below allows you to play with the Ising model (and some variants) by changing the temperature. 

Ernst Ising (actually pronounced Ezing but universally mispronounced as Izing) solved the onedimensional version of his model in his PhD. thesis and didn’t find what he was looking for, namely a phase transition as you see in the watersteam system. He then wrongly inferred that there was no phase transition in two or higher dimensions. Fortunately the model was interesting and simple enough that people continued to investigate it and it was shown that in two dimensions there is indeed a phase transition, though not of the watersteam type. The transition is “weaker”, second order instead of first. In the watersteam system the density suddenly increases as the water condenses at the critical temperature, ie. the boiling point, jumping from almost zero to 1000 kg/m${}^3$. In a second order, Isinglike, phase transition the corresponding parameter, actually the magnetization, goes from zero to a nonzero value smoothly, without a sudden jump.
The model itself is simple. We have a grid, on each grid point lives a spin. Each spin interacts with its neighbours, north, south, east and west. The spins point in only two directions, if two neighbours point in the same direction, $\uparrow \uparrow $ or $\downarrow \downarrow$, the energy is reduced by one, if they point in opposite directions, $\uparrow \downarrow$ or $\downarrow \uparrow$, the energy is increased by one. Slightly more formally, the spins are represented by the variables $s_i$ where $i$ denotes the lattice site. The variables $s_i$ take the values $\pm 1$. The energy is given by, Why does anything interesting happen? Why don’t the spins just all pick a direction and line up, minimizing the total energy? Boltzmann said that the probability of a configuration with energy $E$ is proportional to $g(E)e^{\beta E}$ where $g(E)$ is the number of configurations with energy $E$ and $\beta$ is $1/kT$, where $T$ is temperature and $k$ is a constant with the right units to make the argument of the exponential dimensionless, we will set it to one. At high temperature, small $\beta$, the exponential is approximately $1$, then the state with the largest $g(E)$ is most probable. There are only two states with lowest energy (all $+1$ or all $1$) and it turns for this model there are many, many more configurations with $E = E_{\text{max}}$. At low temperatures, beta is large and the exponential dominates, only by making $E$ small is the probability for a configuration to occur appreciable. We have a tendency to order fighting a tendency to disorder. Energy versus entropy. The Ising model lends itself well to computer simulations. We try to mimic this orderdisorder conflict with the Metropolis algorithm. This says, pick a spin and flip it. If this reduces the energy keep the spin flipped, if it increases the energy keep the spin flipped with probability proportional to $e^{\beta \Delta E}$ , where $\Delta E$ is the change in energy from flipping the spin. Otherwise flip it back. The probabilistic step mimics thermal fluctuations in real systems; the hotter it is the more likely we are to accept flips until whether a spin is up or down is basically random. The colder it is the less likely we are to accept flips that increase the energy and we proceed steadily down to the minimum energy. The most interesting point is the critical point, where order and disorder, energy and entropy, exactly balance. This separates the two (boring) phases, the ‘frozen’ low temperature state and the ‘gaseous’ high energy state. It is a mixture of both, containing bubbles of disorder within order within disorder at infinitum. Play with the animation to see it for yourself. Up spins are red pixels and down spins are green. The $\beta$ slider changes the temperature, the button unit sets all spins equal to $1$, random sets all spins to $\pm 1$ randomly and critical sets the temperature to the critical temperature. The Ising model is the simplest type of spin model and has several generalizations. You can make different shaped lattices; triangular, honeycomb, unionjack. You can make the couplings in the horizontal and vertical directions different, or you can make every coupling on every link different; this type of model is called a spin glass. You can increase the number of dimensions; one and two are exactly solved, three and four are interesting and unsolved analytically; five and higher are somewhat less interesting. This is because of the large number of neighbours each spin has in $N$ dimensions – $2N$. The more neighbours a spin has the more accurate it is to describe the spin as interacting with a fixed background and mean field theory becomes more and more accurate. You could even put several Ising models on the same lattice and make them interact – this is the AshkinTeller model. There are many more. Here I will consider a simple generalization – more than two states, called the Potts model. A spin can point in $N_c$ directions. If two neighbouring spins point in the same direction this decreases the energy, if they point in opposite directions the energy increases, the Ising model has $N_c = 2$. Let’s look at the possible states two neighbours can be in. For the Ising model we have 