Hi!

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:

```
struct sho_functor__ {
std::vector<double>
operator()(const double& t,
const std::vector<double>& y,
const std::vector<double>& theta,
const std::vector<double>& x_r,
const std::vector<int>& x_i, std::ostream* pstream__ = &Rcpp::Rcout) const {
return sho(t, y, theta, x_r, x_i);
}
};
```

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.

Sebastian