Exact solution of the 2D Ising model in an external magnetic field? The 2D Ising model is a thoroughly studied model. One of the remarkable features of the model is that it predicts a hysteresis. However, I cannot seem to find the appropriate literature on this subject. I did do searches on Scopus to find relevant articles. I also found the book "The 2D Ising model" by Barry M. McCoy and Tai Tsun Wu. This is a large book dedicated to the 2D Ising model and has a few paragraphs on hysteresis. One of these shows the expected hysteresis: a "smooth" flip from the -1 to 1 normalized magnetization state.
However, this is only for the first row of the lattice. This is not what I am searching. Is there an analytical expression for the whole net magnetization of the whole lattice, in an external magnetic field? The field is swept from some negative value to some positive value, and then from positive back to negative, while keeping the temperature constant.
A pointer to a numerical solution or to a simulation is fine too.
EDIT:
Background
I am using the Metropolis algorithm to simulate an Ising lattice. See the following figure generated by the model:

Apologies for the weird legend. What you see: the model starts at T=0, then is swept in linear steps to T=2 in 100 steps. Each step sweeps 100 times over the lattice. The lattice is 20x20 and uses periodic boundary conditions. I'm using J = 1 and k = 1, where J is the coupling constant between the spins. This first step is done to generate an appropriate "initial condition" before the field is turned on.
Then the field is swept from -0.1 to 1, in 1000 steps, linearly. The same is done but now backwards.
As you can see, there is an hysteresis as expected. However, I would expect the hysteresis to be "smooth", but it appears to be a step function (or "square" hysteresis).
I tried zooming in on the hysteresis region, but it either moves or gets even sharper (the intermediate points, which appear to be artifacts, are deleted).
I'm therefore wondering wether this hysteresis is what I should see, or that I'm doing something wrong in my simulation. Maybe this outcome is a wrong prediction when using the Metropolis algorithm?
 A: The Ising model is indeed very interesting!
In 2-dimensions, there is an analytical solution, in the case of no applied field.  It is very complicated, and when it first came out, it consisted of 30 pages of very challenging maths.  Only to be 'simplified' down to about 15 pages in the 60s.
With an applied field, or in higher dimensions, it is typical to numerically simulate the Ising model.  The system is very amicable to Monte Carlo methods.  You can look up the Metropolis algorithm.  I found a very useful resource was the book Monte Carlo Methods in Classical Statistical Physics by Janke.   Here is an excerpt https://www.physik.uni-leipzig.de/~janke/Paper/lviv-ising-lecture-janke.pdf
A lot of very interesting physics can be seen in these numerical simulations.  Here's one of many online simulations of the 2D Ising Model:
http://physics.weber.edu/schroeder/software/demos/IsingModel.html
If you want an offline simulation, you can find some (along with many other interesting simulations) at NetLogo http://ccl.northwestern.edu/netlogo/
Edit:
Your hysteresis curve is definitely believable, if you are under $T_\textrm{C}$, which you easily are.  One of the results from Onsager's original paper, is that the spontaneous magnetisation, while below $T_\text{C}$, is given by:
$$|m|=\Bigg(1-\sinh\Big(\frac{2J}{T}\Big)^{-4}\Bigg)^\frac{1}{8}.$$
I've included an image of when I was simulating the 2D ising model, and in the inset, the red line is the result from the exact solution.  (No applied field.)
                                  

What you can see from this, is that when you are below the critical temperature, you can expect a large spontaneous magnetism.  At $T=1$, the system in equilibrium is nearly fully magnetised, spontaneously.  Although there isn't an exact solution, you can imagine what will happen when you apply an external field.  If the field is in the direction of your magnetisation, it will slightly increase the magnetisation.  (If the field is strong enough, it will become fully magnetised.)  On the other hand, if the field is small, and in the opposite direction, it will, after a quick period, entirely flip the spins, so they are at a little above their spontaneous magnetisation value, in the opposite direction.
I've mocked up a hysteresis curve based on this reasoning, at $T=1$ (note the heavily discontinuous scale!).  The spontaneous magnetisation value of $m=0.999$ comes from the formula I quoted earlier.  At $H=0$, both $m=\pm0.999$ are equilibrium states, but as soon as you slightly alter $H=\pm\delta$, where $\delta$ is small, you pick one of them.
                                  
There's one crucial thing missing from this logic, and that's the problem of nucleation energy costs.  If all your spins point one way, the system can exist out of equilibrium, due to it being hard to flip all the spins.  For all the spins to flip, requires at first a small cluster of spins to flip against their neighborhood, which is energetically unfavourable, despite the flip being with the field.  That leads to a curve with hysteresis behaviour, and it's in good agreement with what you've seen.  Once a field is reached strong enough to overcome the nucleation barrier, the entire grid of spins can flip.
                                  
A: There exist analytic results for the second derivative of the of the free energy w.r.t. the magnetic field at field strength zero (i.e. the susceptibility), see e.g. here. These results were obtained using experimental mathematical work where one guesses exact expressions based on series expansions. There is no exact expression of the full susceptibility, but it can be expressed as an infinite sum of D-finite functions (functions that satisfy linear differential equations with polynomial coefficients). 
