I've been using the odeint library for some time now - and I really like it, thanks for the great work! I know that the odeint library is currently being reviewed but I'm actually writing here to ask a couple of questions about my usage scenario, hope that I'm in the right place.

My usage scenario involves neural networks where I have a set of n equations for the N neurons themselves and another m equations for all the M synapses (connections) where usually M >> N and the two sets of equations depends on each other. What I'm trying to do is to put all these N + M equations into one large system and to integrate them at once while also keeping the state vectors separated. This is because of the different dimensionality of the two sets of equations and because my intention is to later use OpenCL to do the calculations (which as I understand works better if you use something like a vexcl::multivector<double, M> + another vexcl::multivector<double N> and I don't see how to pack all of this in one large state vector);

So my first question is: can I somehow "combine" two sets of state vectors into one thing or integrate two systems in parallel (synchronously)? I tried the former (one combined state object) and even defined my own combined algebra (a std::pair of std::vectors - later vexcl::multivectors) but then stumbled into problems with my value type (which I defined as a std::array<double, 2>). It seems that I have to define all possible mathematical and comparison operations for my pairs (arrays).

And my second question: my problem involves integrating the system until one of the neurons fires - at which point it's state variable (membrane potential) jumps to a fixed value, i.e., I have a sharp discontinuity in the state variable. My way of dealing with this right now was to implement a custom integration function (based on the adaptive integration one) which has a built in function that checks if the threshold was reached at every integration step. When this happens I break the integration, reset the state of the neuron in question to the fixed value and start integrating again. I wonder if it's possible to solve this in an easier way... Is it possible for example to change the state on the fly and can the system handle discontinuities? (i.e. can I integrate without stopping each time a neuron fires).