Exact solution of the 2D Ising model in an external magnetic field?

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.