Ok, so here is what I think is really killing the performance: The problem as written uses a piece-wise linear approximation to the infections function. For 63 time-intervals a slope and an offset is being fitted. So we have 2x63 parameters just for that… and these are all passed into the ODE RHS. The problem is that per parameter we get again N ODE’s - here N is the number of states, which is 8 in this case… so we end up having to solve just for that about 900 ODEs!!!
BUT… we can completely avoid this issue by doing the integration in steps along the 63 time-intervals. So integrating by time-interval 63x we can simple pass for each ODE solve one offset and one slope…that reduces the extra number of ODEs per solve from 900 to just 18. It’s a bit of a hassle to get the indexing right (and I messed it up in my first attempt), but hopefully I find the time to come up with a working model.
The bottom line here is: NEVER pass in parameters into the ODE RHS unless you have to and really work hard to avoid it.
Ah… and adjoint ODE will be interesting here to try, of course.