What is an integrable system? What is an integrable system, and what is the significance of such systems? (Maybe it is easier to explain what a non-integrable system is.) In particular, is there a dichotomy between "integrable" and "chaotic"? (There is an interesting Wikipedia article, but I don't find it completely satisfying.)
Update (Dec 2010): Thanks for the many excellent answers. I came across another quote from Nigel Hitchin:
"Integrability of a system of differential equations should manifest
itself through some generally recognizable features:

*

*the existence of many conserved quantities


*the presence of algebraic geometry


*the ability to give explicit solutions.
These guidelines would be interpreted in a very broad sense."
(If there are some aspects mentioned by Hitchin not addressed by the current answers, additions are welcome...)
Closely related questions:

*

*What does it mean for a differential equation "to be integrable"?

*basic questions on quantum integrable systems
 A: Since, it hasn't been mentioned yet a short addition to José Figueroa-O'Farrill's answer. I will only talk about the finite dimensional case. So let's assume that $\dim(P) = 2n$. Then the Hamiltonian flow is integrable if there exist these $n$ functions $f_1, \dots, f_n$ which are in involution with respect to the Poisson structure.
Now, the cool thing is that there exist action angle coordinates. These means we can conjugate our possibly complicated dynamics to the simple dynamics
$$
 \partial_t I_j = 0,\quad \partial_t \theta_j  = I_j,\quad j=1,\dots,n 
$$
this is something, we can all solve since it is just linear. Note: We will have $I_j  = f_j(\text{orbit})$, which is time independent.
As a possible application, KAM theory is usually formulated as an application to systems in action angle coordinates. This in turn implies that integrable systems are stable (in a subtle measure theoretic sense). But I think this is what is meant with "integrable $\neq$ chaos". We have a great form of perturbation theory for integrable systems.
A: After reading several books and articles about integrable systems, and after several years of work in the field, I consider particularly meaningful  the following quotation from  Frederic Helein's book 'Constant mean curvature surfaces, harmonic maps and integrable systems', Lectures in Mathematics, ETH Zurich, Birkhauser Basel (2001):
"...working on completely integrable systems is based on a contemplation of some very exceptional equations which hide a Platonic structure: although these equations do not look trivial a priori, we shall discover that they are elementary, once we understand how they are encoded in the language of symplectic geometry, Lie groups and algebraic geometry. It will turn out that this contemplation is fruitful and lead to many results"
A: I'll take off from the questioner's suggesting that maybe it's better to say what is a NON-integrable system is.
The Newtonian planar three body problem, for most masses, has been proven to be non-integrable.
Before Poincare, there seemed to be a kind of general hope in the air that every autonomous Hamiltonian system was integrable.
One of Poincare's big claims to fame, proved within his Les Methodes Nouvelles de Mecanique Celeste, was that the planar three-body problem is not completely integrable.   It is the  dynamical systems equivalent to Galois' work on  quintics. 
Specifically, Poincare proved that besides the energy, angular momentum and linear momentum there are no other ANALYTIC functions on phase space
which Poisson commute with the energy.  (To be more careful:
any 'other' such function is a function of energy, angular momentum, and linear momentum.
And his proof, or its extensions, only holds in the parameter region where
  one of the mass dominates the other two.  It is still possible that for very
  special masses and angular momenta/ energies the system is integrable.
  No one believes this.)  As best I can tell, existence of additional   smooth
  integrals (with fractal-like level sets) is still open, at least in most cases. 
Poincare's impossibitly proof is based on his  discovery of   what is nowadays called a "homoclinic tangle"
  embedded within the restricted three body problem,  viewed in a rotating frame.
  In this tangle, the unstable and  stable manifold of some point  (an orbit in the non-rotating
  inertial frame)    cross each other
   infinitely often, these crossing points having the point in its closure.
  Roughly speaking, an additional integral would have to be constant  along this complicated set.
Now use  the fact that if  the zeros of an analytic function have an accumulation point then that  function is zero
to conclude that the function is zero. 
Before Poincare (and I suppose since) mathematicians and in particular astronomers   spent much energy searching for  sequences of changes of variables which made the system "more and more integrable".  Poincare realized the  series defining their transformations were divergent -- hence his interest in divergent series.
This divergence problem is the "small denominators problem" and getting around it by putting number theoretic conditions on frequencies appearing  is at the heart of the KAM theorem.   
A: This is, of course, a very good question.  I should preface with the disclaimer that despite having worked on some aspects of integrability, I do not consider myself an expert.  However I have thought about this question on and (mostly) off.
I will restrict myself to integrability in classical (i.e., hamiltonian) mechanics, since quantum integrability has to my mind a very different flavour.
The standard definition, which you can find in the wikipedia article you linked to, is that of Liouville.  Given a Poisson manifold $P$ parametrising the states of a mechanical system, a hamiltonian function $H \in C^\infty(P)$ defines a vector field $\lbrace H,-\rbrace$, whose flows are the classical trajectories of the system.  A function $f \in C^\infty(P)$ which Poisson-commutes with $H$ is constant along the classical trajectories and hence is called a conserved quantity.  The Jacobi identity for the Poisson bracket says that if $f,g \in C^\infty(P)$ are conserved quantities so is their Poisson bracket $\lbrace f,g\rbrace$.  Two conserved quantities are said to be in involution if they Poisson-commute.  The system is said to be classically integrable if it admits "as many as possible" independent conserved quantities $f_1,f_2,\dots$ in involution.  Independence means that the set of points of $P$ where their derivatives $df_1,df_2,\dots$ are linearly independent is dense.
I'm being purposefully vague above.  If $P$ is a finite-dimensional and symplectic, hence of even dimension $2n$, then "as many as possible" means $n$.  (One can include $H$ among the conserved quantities.)  However there are interesting infinite-dimensional examples (e.g., KdV hierarchy and its cousins) where $P$ is only Poisson and "as many as possible" means in practice an infinite number of conserved quantities.  Also it is not strictly necessary for the conserved quantities to be in involution, but one can allow the Lie subalgebra of $C^\infty(P)$ they span to be solvable but nonabelian.
Now the reason that integrability seems to be such a slippery notion is that one can argue that "locally" any reasonable hamiltonian system is integrable in this sense.  The hallmark of integrability, according to the practitioners anyway, seems to be coordinate-dependent.  I mean this in the sense that $P$ is not usually given abstractly as a manifold, but comes with a given coordinate chart.  Integrability then requires the conserved quantities to be written as local expressions (e.g., differential polynomials,...) of the given coordinates.
A: I have found an article "Quantum signatures of an integrable system with a chaotic scattering map" here:
http://www.iop.org/EJ/abstract/0305-4470/28/6/008
So apparently some integrable systems can have chaotic scattering maps.
A: For a Hamiltonian system, integrable means the solution lies on a surface (in phase space).  The dimension of the surface depends on how many integrals there are.  And as long as we are talking about Hamiltonian systems, chaotic and integrable are indeed complements of each other, but this is not the case for dynamical systems in general. 
A: The simple answer is that a $2n$-dimensional Hamiltonian system of ODE is integrable if it has $n$ (functionally) independent constants of the motion that are "in involution". (Functionally independent means none of them can be written as a function of the others. And "in involution" means that their Poisson Brackets all vanish -- a somewhat technical condition I won't define carefully (* but see below), but instead refer you to: http://en.wikipedia.org/wiki/Poisson_bracket). The simplest and the motivating example is the $n$-dimensional Harmonic Oscillator. What makes integrable systems remarkable and interesting is that one can find so-called "action angle variables" for them, in terms of which the time-evolution of any orbit becomes transparent. 
For a more detailed and modern discussion you may find an expository article I wrote in the Bulletin of The AMS useful. It is called "On the Symmetries of Solitons", and you can download it as pdf here:
http://www.ams.org/journals/bull/1997-34-04/S0273-0979-97-00732-5/
It is primarily about the infinite dimensional theory of integrable systems, like SGE (the Sine-Gordon Equation), KdV (Korteweg deVries) , and NLS (non-linear Schrodinger equation), but it starts out with an exposition of the classic finite dimensional theory.


*

*Here is a little bit about what the Poisson bracket of two functions is that explains its meaning and why two functions with vanishing Poisson bracket are said to "Poisson commute". Recall that in Hamiltonian mechanics there is a natural non degenerate two-form $\omega = \sum_i dp_i \wedge dq_i $. This defines (by contraction with $\omega$) a bijective correspondence between vector fields and differential 1-forms. OK then -- given two functions $f$ and $g$, let $F$ and $G$ be the vector fields corresponding to the 1-forms $df$ and $dg$. Then the Poisson bracket of $f$ and $g$ is the function $h$ such that $dh$ corresponds to the vector field $[F,G]$, the usual commutator bracket of the vector fields $F$ and $G$.
Thus two functions Poisson commute iff the vector fields corresponding to their differentials commute, i.e., iff the flows defined by these vector fields commute. So if a Hamiltonian vector field (on a compact $2n$-dimensional symplectic manifold $M$) is integrable, then it belongs to an $n$-dimensional family of commuting vector fields that generate a torus action on $M$. And this is where the action-angle variables come from: the level surfaces of the action variables are the torus orbits and the angle variables are the angles coordinates for the $n$ circles whose product gives a torus orbit.

A: "All integrable Hamiltonian systems are alike, while each nonintegrable one is nonintegrable in its own way", Valerij V. Kozlov (after L. N. Tolstoi of course), Symmetries, Topology and Resonances in Hamiltonian Mechanics, Springer, 1996
A: I don't think that one could  say that there is a dichotomy between integrable and chaotic systems. There is certainly a huge chunk in the middle. By a chaotic system we often mean a system where trajectories of points deviate exponentially with time, a canonical example is the Arnold (or Anosov) cat's map. In this case a generic trajectory is of course everywhere dense in the phase space. This is related to ergodicity (in the case when there is a measure preserved by the system). But of course not every ergodic system is chaotic. There are different degrees of chaos, mixing, strong mixing, etc.
On the contrary for an integrable system the motion of every trajectory is quasi-periodic, it stays forever on a half-dimensional torus, such systems are rare. A little perturbation of such a system is not integrable anymore. KAM theory describes the residue of integrability of the perturbation, while Arnol'd diffusion is about trajectories that don't move quasiperiodically anymore. 
There is one amazing example due to Moser, that shows how the cat's map can "happen"  on a degenerate level of an integrable system: page 6 in
http://arxiv.org/PS_cache/arxiv/pdf/0810/0810.5713v1.pdf
A: The above answers deal mostly with finite-dimensional systems. As for the (systems of) PDEs, you typically need the Lax pair or a zero curvature representation (see e.g. the Takhtajan--Faddeev book mentioned in the wikipedia entry you linked to for the definition of the latter) or something else like that. To the best of my knowledge, the complete understanding of what is an integrable system for the case of three (3D) or more independent variables is still missing. In particular, for the case of three independent variables (a.k.a. 3D or (2+1)D) the overwhelming majority of examples are generalizations of the systems with two independent variables. These generalizations are constructed using the so-called central extension procedure (e.g. the KP equation is related to KdV in this way). Many integrable partial differential systems in three independent variables and apparently the overwhelming majority thereof in four or more independent variables are dispersionless a.k.a. hydrodynamic-type, i.e., can be written as first-order homogeneous quasilinear systems, see e.g. this paper and references therein for details. 
As for the reading suggestions, in addition to the Takhtajan--Faddeev book cited above, you can look e.g. into a fairly recent book Introduction to classical integrable systems by Babelon, Bernard and Talon, and into the book Multi-Hamiltonian theory of dynamical systems by Maciej Blaszak which covers the central extension stuff in a pretty straightforward fashion. Both books have extensive bibliographies with further references to look into. 
Now, as for classification and identification of (new) integrable systems of PDEs, at least in two independent variables, it turns out that the (infinitesimal higher) symmetries play an important role here.
A recent collective monograph Integrability, edited by A.V. Mikhailov and published by Springer in 2009, could be a good starting point in this direction. See also another recent book Algebraic theory of differential equations edited by MacCallum and Mikhailov and published by Cambridge University Press. 
For a general introduction to the subject of symmetries of (systems of) PDEs, I can recommend the book Applications of Lie groups to differential equations by Peter Olver.
A: The integrability conditions for the existence of a Lagrangian or Hamiltonian are known as "conditions of variational self-adjointness."
The conditions are studied in the context of the Inverse Problem, which


*

*Ruggero Maria Santilli, Foundations of Theoretical Mechanics: The
Inverse Problem in Newtonian Mechanics, vol. 1, 2 vols. (New York:
Springer-Verlag, 1978), p. 10,


formulates as:

Given the totality of solutions $y(x) = \left\{y^1(x), \ldots, y^n(x)\right\}$ of a system of $n$ ordinary differential equations of order $r$,
  $$F_k\left(x, y^{(0)}, y^{(1)}, \ldots, y^{(r)}\right)=0\qquad\text{(I.23)}$$
  $$y^{(i)}=\frac{d^iy}{dx^i},\qquad i=1,\ldots,r,\qquad k=1,2,\ldots,n,$$
  determine whether there exists a functional
  $$A(y)=\int_{x_1}^{x_2}dxL\left(x,y^{(0)},\ldots,y^{(r-1)}\right)\qquad\text{(I.24)}$$
  which admits such solutions as extremals.

It appears Helmholtz was the first to study the Inverse Problem (ibid., p. 12):

The necessary and sufficient conditions for the existence of a solution $L$ of system (I.29)* were apparently formulated for the first time by Helmholtz
  (1887)26 on quite remarkable intuitional grounds. In essence, Helmholtz's starting point was the property of the self-adjointness of Lagrange's equations, i.e., their system of variational forms coincides with the adjoint system (see Chapter 2 and following). This is a property which goes back to Jacobi (1837).27 Without providing a rigorous proof, Helmholtz indicated that the necessary and sufficient condition for the existence of a solution $L$ of system (I.28)** is that the system $F_k = 0$ be self-adjoint.
  26. Helmholtz did not consider an explicit dependence of the equations of motion on time. Subsequent studies indicated that his findings were insensitive to such a dependence.27.  The equations of variations of Lagrange's equations or, equivalently, of Euler's equations of a variational problem, are called Jacobi's equations in the current literature of the calculus of variations. We shall use the same terminology for our Newtonian analysis.

*(I.29) is the Euler-Lagrange equation corresponding to (I.24) when $n>1$, $r=2$.**(I.28) is the case when $n=r=1$.
This same analysis of conditions for variational self-adjointness of a Lagrangian can be applied to Hamiltonians, as Hamiltonians are simply the Legendre transform of Lagrangians (cf. Callen's Thermodynamics and an Introduction to Thermostatistics §5.2 (pp. 137-145) for a good introduction to Legendre transforms).
Helmholtz (1887) is:


*

*Ueber die physikalische Bedeutung des Prinzips der kleinsten Wirkung ("On the physical meaning of the principle of least action")J. Reine Angew. Math., 100 (1887), pp. 137–166


*

*See esp. §3 "The derivation of the kinetic potential from the value of energy." (pp. 18-21 of the English translation), where he (as he says in his outline on p. 5) "treats the opposite problem, namely, that of deriving $H$ [the Hamiltonian] from $E$ [the energy]."


A: This is soft -- but I think of an integrable system as one whose dynamics are dominated by algebra.  For finite dimensional integrable systems, the symmetries (related to conserved quantities by Noether's theorem) force the trajectories to live on half-dimensional tori.  For infinite dimensional integrable systems, where the flow on the scattering data is isospectral the symmetries force solutions to be n-soliton solutions plus dispersive modes.  
There is a blog post of Terry Tao's (apologies for not having the link) which talks about how algebra is the right tool to understand structure while analysis is the right tool to understand randomness.  The claim is that one mark of an good problem is the presence of an interesting relationship between structure and randomness and hence the requirement that both algebra and analysis be used -- to some degree -- in order to get a good answer to the problem.  The soliton resolution conjecture is by this standard a good problem because the asymptotic n-soliton solutions are fundamentally algebraic while the dispersive modes are fundamentally analytic objects.
I agree with Dmitri that there isn't a dichotomy.  The symmetries can have a large or small role in the dynamics as can the ergodicity.
A: I'll give a bit of a physics definition.  (Reference is "A Brief Introduction to Classical, Statistical and Quantum Mechanics" by B\"uhler.)
"A mechanical system is called integrable if we can reduce its solution to a sequence of quadratures."
So, literally, an integrable system (in this view) is one that can be solved by a sequence of integrals (which may not be explicitly solvable in elementary functions, of course).  To connect to other answers, this should only work out when there are enough symmetries for us to write down and integrate.
A: I would like to add one more example of integrability which refers to Hopf algebras and is probably the easiest to formulate. It naturally arises in spin chain physics, but can be treated abstractly as well. Consider (semi-simple) Lie algebra $\mathfrak{g}$, its universal enveloping algebra $\mathfrak{h}=U(\mathfrak{g})$ and its Hopf double. The latter has coproduct homomorphism $\Delta: \mathfrak{g}\to \mathfrak{g}\otimes\mathfrak{g}$. Now let us consider an operator $\mathfrak{R}$ (so-called, R-matrix) as the following mapping $\mathfrak{R}:\mathfrak{h}\otimes\mathfrak{h}\to\mathfrak{h}\otimes\mathfrak{h}$, meaning that $\mathfrak{R}$ is some tensor product of polynomials of elements from $\mathfrak{g}$. The integrability condition reads
$[\Delta,\mathfrak{R}]=0,$
viz. the coproduct should commute with the R-matrix. It is now a matter of several lines of simple calculations to show that the Yang-Baxter equation on the R-matrix, which is frequently referred to as the necessary condition for integrability follows [see, e.g. Kassel "Quantum Gorups"]. 
