I think everything that has been discussed in the CoDatMo thread should be relevant again.
No matter which solver you use, as soon as you have non-constant (over time) parameters, it is very very (very) likely that you have to split up the solution of the ODE by hand into separate time intervals in which the ODE parameters are smooth. And don’t do linear or binary or whatever search to identify the correct beta[i].
Also, don’t forget about finding an efficient representation of the effects you want to capture.