The Ising model is a model, originally developed to describe ferromagnetism, but subsequently extended to more problems.
Basically, it is an interaction model for spins. Imagine you have a system which is a collection of $N$ spins. Each spin $S_i$ has two possible states $+1$ or $-1$. Here you can imagine already a possible extension to more states. You can also imagine a different interpretation to the spins: $-1$ is a box containing no gas particle, $+1$ is a box containing a gas particle.
But I'm getting ahead of myself. Let's stick to spins for now.
The next step is to define the energy of the system.
$$E= - \sum_i h_i S_i - \sum_{i \neq j} J_{ij} S_i S_j$$
The first term can be interpreted as the contribution to the energy of the interaction of a spin with a local magnetic field. If the magnetic field is the same for all spins, then $h_i=h$ for all $i$. You see that aligning with the field will means a lower energy than going against the field.
The second term represents interactions between spins within the system. If $J_{ij}>0$ spins which align will contribute negatively to the energy of the system, thus lowering the total energy. If $J_{ij}<0$ then anti-alignment will contribute.
What I still haven't specified is how the spins are structured. By choosing the coefficients $J_{ij}$ appropriately, I can introduce this structure. Suppose I want a 1-dimensional system, which is the one Ising originally solved. Then, you have a countable infinity of spins arranged along the real line with equal spacing. Ising imposed an interaction only allong neighboring spins. So spin $S_n$ can interact with spin $S_{n-1}$ and spin $S_{n+1}$. The energy formula I gave above becomes:
$$E= - J \sum_{n} S_n S_{n+1} \; ,$$
if all spins can interact equally strongly.
Now, to each configuration of the system corresponds a certain energy. In statistical mechanics, we know that the probability of a certain configuration is
$$P(\{S\}) \sim e^{-E(\{S\})/kT}$$
where $T$ is the temperature. Or, if we compute the partition sum
$$Z=\sum_{\{S\}} e^{-E(\{S\})/kT}$$
we can deduce the complete equilibrium thermodynamic properties of the system. In particular, we can see if there are phase transitions for certain values of the parameters. (There is none in the 1D case, which at the time, combined with the lack of attention his model was getting, made Ising abandon physics.)
Of course, there are many more generalizations of this model. You can also make dynamical versions of the model where you are not just interested in the equilibrium configurations.
Here's the wikipedia page for Ising models.
Some references:
H.E. Stanley 'Introduction to Phase Transitions and Critical Phenomena' Clarendon Press Oxford
J.M. Yeomans 'Statistical Mechanics of Phase Transitions' Clarendon
Press Oxford
R.J. Baxter 'Exactly Solved Models in Statistical Mechanics' Academic New York