I discovered that functions which have calls to integrate_ode_bdf cannot be exported to R for a very subtle reason: the expose_stan_functions functions throws away all the template magic wrt to var. Therefore, the ODE functor gets turned into something like:

Therefore, the bdf solver cannot anymore obtain any Jacobian wrt to the states which are required for the stiff solver. Hence, the compilation aborts with a nasty compilation error. Anything we can do about this?

If the answer is yes, I am happy to file an issue to rstan, but I am not quite sure.

Hmmâ€¦ templates are not really needed. What we would need is an additional functor which has the states as vector of vars instead of doubles (or insider the functor a function with this signature). I suppose this is not really practical.

Any suggestions what we should do:

Nothing (probably I am the only person noticing)

via a regex replace integrate_ode_bdf with integrate_ode_rk45 and warn?

file an issue for the record to fix at some point when we have better ideas?

I am not totally sure what the best option would be, but have you evaluated whether including a little bit of C++ would solve your problem for the time being?

Yes, including some code would work. That code is not trivial to get, but doable for me (supplying the Jacobian in analytic form). Hence, this solves the problem for me, but not for general users.