Hi Brian, I know we've talked about this before, but here are some more general suggestions about your problem.
So, I changed the implementation so that everything was done in "one-shot" (luckily the only non-DE rule I had so far was easily converted to a DE rule, although there is one coming up which I don't see any easy conversion for...[that is the subject of this question](http://scicomp.stackexchange.com/questions/14765/scipy-integrate-odeint-how-...)).
Having looked at that, your problem is not so bad. You will need to solve a piecewise (a.k.a. "hybrid") system, but you shouldn't be discretizing it based on time but based on events -- i.e. when your particle crosses into a new discrete area you are calling p. You then have a discrete state update for your parameters and then you can restart your smooth integrator until the next transition. This is a well-studied problem in this form and is numerically soluble with a hybrid solver such as that provided by PyDSTool. Between positional transitions, your dynamical system is smooth and can be solved with a regular ODE solver. The Jacobian within each domain should be simple enough to pre-calculate (as a function of p or <p>). There are a few hybrid model examples with PyDSTool on the tutorial and I am willing to help a bit with the setup once you've given a new script for your problem a shot. Take a copy of a helpful example (e.g. the SLIP pogo stick dynamics) and adapt it to set up what you can. Put in some comments etc. of what you need to happen. PyDSTool will be able to convert the ODEs to C code automatically but not the transition rules, which will still happen in python. This will not be the *fastest* way to solve it but, more importantly, this way will give you an accurate solution in a form that you can understand and manipulate, and you can worry about optimizing speed later, IMO. -Rob