Stage-wise controls for IRK methods

A big thanks to @vmerc !
- stagewise control is now the default for GL2 and GL3 scheme, with the stepwise control still available through options `gauss_legendre_2|3_constant_control`
- [DONE] manual sparsity patterns are not available yet, but should not be too difficult to adapt from the stepwise control formulation
- solution currently outputs the 'averaged' control at each step, not the individual stage controls (which could be retrieved with a bit of work in the high-level generic getter). The averaged control uses the same weights `b_j` as the slopes in the discretized dynamics `x_i+1 = x_i + h_i sum_j b_j f_ij`. Recall that the `f_ij` are part of the NLP unknown, and subject to the equations `f_ij = f(t_ij, x_ij, u_ij)`. The state values at the stages are treated as internal variables and computed via `x_ij = x_i + sum_k a_jk f_ik`. Alternately, it may be possible to use the `x_ij` as NLP unknowns and compute the `f_ij` instead with proper buffering. To be investigated.
- another related question would be the handling of path constraints, which are currently evaluated at time steps using the first stage controls, ie `g(t_i, x_i, u_i1)`. We could use the averaged control instead of the first one, or conversely evaluate the constraint at the stage times, reusing the `x_ij` from the dynamics part.

@ocots @jbcaillau @joseph-gergaud this version will be the 1.0.9 release