I must say I am like Colin and tend to think of the two phenomenons as overlapping and the way I would describe diffraction would be roughly as Colin does (which maybe isn't too surprising given that our backgrounds overlap). I must admit to never really having cared about the difference deeply - but now I read @dmckee 's definition, I think I like it because it is cleanly discriminating between the two and moreover when I think about it it does pretty much describe my own gut definition if I had to give one.
I think I can reword @dmckee 's diffraction definition as: "how a propagating light field, having encoded perturbations imposed on it at one transverse plane (i.e. orthogonal to the nominal propagation direction), and expresses these perturbations at another transverse plane".
Here's my own understanding. I think of diffraction mostly as:
A wavefront scrambling process arising from the divergence of a light field's constituent plane waves
Consider a field on a plane, say $z = 0$ and split it up using Fourier decomposition of the field variation over the plane $z=0$ into constituent plane waves, which are "modes" of Maxwell's equations insofar that their propagation description is simply that the fields become phase delayed by a simple scale factor $\exp(i\,\mathbf{k}\cdot\mathbf{\Delta r})$ under the action of a translation $\mathbf{\Delta r}$. Each constituent plane wave has a different direction defined by the wavevector $\left(k_x, k_y, k_z\right)$ with $k^2 = k_x^2 + k_y^2 + k_z^2$ (i.e. the Fourier space equivalent of Helmholtz's equation), that is, all the wavevectors have the same magnitude but different directions. So, when we ask what the field looks like at a different value of $z$, we build the field up from our plane wave constituents at this point (use an inverse Fourier transform). However, now, because the wavevectors are all in different directions, the plane waves have all undergone different phase delays in reaching the new value of $z$ (even though their phase advances by $k$ radians per unit length in the direction of the respective wave vector). Therefore, the field's configuration gets scrambled by all these different phase delays. I sketch this idea in a drawing below:
To illustrate this idea further, we think of a one-dimensional problem, so we have a uniformly lit slit of some finite width $w$ modeling the laser output; in this simplified system that there are only 2D wave vectors. The screen with the slit is in the $z = 0$ plane and the one orthogonal direction is the $x$ axis. All the Cartesian components of the fields fulfil the same (Helmholtz) equation, so we can discuss the principles by just looking at one scalar field $\psi$ (say, the electric field's $x$-component). Each plane wave has the form $\psi(k_x) = \exp\left(i \,(k_x\, x + k_z\, z)\right)$ The Fourier transform of the field output from the slit is then (I'll leave out factors of $2\pi$ in the unitary FT because scale factors don't affect the following):
$$\frac{\sin\left(\frac{w\, k_x}{2}\right)}{k_x} \quad\quad\quad(1)$$
where $w$ is the slit width, and unless the slit is very wide, the Fourier transform has a wide spread of frequencies. This means that for $z = 0^+$ ("immediately downstream" of the slit's output) the field is the superposition
$$\int\limits_{-\infty}^\infty \frac{\sin\left(\frac{w\, k_x}{2}\right)}{k_x} \exp\left(i\, (k_x\, x + k_z\, z)\right) \mathrm{d} k_x\quad\quad\quad(2)$$
When we plug $z = 0$ in, the integral is simply the inverse FT of (1) and we get our original slit field. But now put some nonzero value of $z$ in: because $k_x^2 + k_z^2 = k^2$, we have $k_z = \sqrt{k^2 - k_x^2}$ (assuming the field is running in the $+z$ direction), we get
$$\int\limits_{-\infty}^\infty \frac{\sin\left(\frac{w\, k_x}{2}\right)}{k_x} \exp\left(i\, (k_x\, x + \sqrt{k^2 - k_x^2}\, z)\right)\, \mathrm{d} k_x\quad\quad\quad(3)$$
You can see the "scrambling", $k_x$-dependent phase factor $\exp(i\, \sqrt{k^2 - k_x^2}\, z) = \exp\left(i\, k\, \cos\theta_x\,\right)$ (where $\theta_x$ is the angle that the plane wave with wavevector $(k_x, k_z)$ makes with the $z$-axis) will yield the complicated scrambling you see as "diffraction".
A phenomenon I would definitely have preferred to call "interference" rather than "diffraction" would be multipathing: i.e. the splitting of a beam of light into two distinct beams e.g. for interferometry and the patterns of "interference" one gets when the beams are brought back together again.
To look at a common "ambiguous" situation: one speaks equally of "single slit interference" or "single slit diffraction" for the same thing: given what I have said above, I would rather call this "diffraction".
Lastly: Ultimately diffraction can be thought of as a special case of interference. Instead of spatially separated beams interfering, however, we take heed that the diffracted wavefront arises from the interference between the separate plane waves that make up a light field and because the pathlength for these interfering "beams" depends on propagation direction, you get the wavefront "scrambling" I spoke of.
