Ah, I see. I think I may just try my luck with ode_ckrk for now then before doing any kind of manual tinkering.
Also, w.r.t. inference of multidimensional ODE states, would it be possible for me to do
y ~ normal(y_hat, sigma)
in the model block, where y and y_hat are multidimensional array[state_dim, N] real, or do I need to melt to a univariate structure? I recall from a couple of years ago that this vectorization was not possible with matrix types, but did not know if the operation was possible for a multidimensional array of reals.
Thanks again!