It is a method for solving ordinary differential equations. Its origins are in something called a LaPlace transform.
The LaPlace transform is defined as
L{ f(t) } = integral (0 to inf) f(t) * exp( -st) dt = F(s).
Note that when you integrate with respect to t with bounds, it is no longer a function of t, but rather a function of s. s has units of inverse of whatever your space is -- if t is a distance, s is an inverse distance; if t is time, s is inverse time (aka frequency).
There are a lot of interesting properties of the LaPlace transform. One in particular is the transform of a derivative is the same as the transform of the original times s minus the function at t = 0, or
L{ df(t)/dt } = s*F(s) - f(0).
This is really interesting because it means if you take the LaPlace transform of a differential equation, it becomes an algebraic expression that can more easily be solved. After solving in terms of s, typically the inverse transform is found by partial fraction decomposition and looking up the individual terms on a table.
When you solve a differential equation in which there is an unknown forcing function, such as
d^2 f(t)/dt^2 + a*df(t)/dt + bf(t) = g(t)
and you take the LaPlace transform of it and solve for F(s), you end up with something like this:
F(s) = [rational function in s] * G(s)
Now how do you handle the inverse LaPlace transform when you have multiplication like that? It turns out [with a proof that isn't that hard] that the inverse LaPlace transform of two multiplied functions is the two functions convolved. This gives you a generalized time-domain (or normal-space domain) solution to your differential equation with unknown forcing function.
It is used very often, given that with a computer convolution can be done in real-time whereas transform based analysis requires a significant delay. In circuit simulation, filters and amplifiers that have non-flat frequency response are described by differential equations, and if you want to look at transient circuit behavior or the response of a circuit to something that isn't a sine wave in time rather than in frequency, it is computed using convolution.
There is also a method that arrives at the use of convolution for solving differential equations without touching the LaPlace of Fourier transforms, but it's more difficult to describe as it requires the use of Dirac Delta functions and explaining LTI systems.