FWIW, here’s full model based on the above ODE and *new* solver interface. The model is adapted from my StanCon UK tutorial https://github.com/metrumresearchgroup/torsten_stancon_2019_tutorial/tree/master/RScripts/model/chemical_reactions

```
functions{
vector reaction(real t, vector x, real[] p){
vector[3] dxdt;
real p1 = p[1];
real p2 = p[2];
real p3 = p[3];
dxdt[1] = -p1*x[1] + p2*x[2]*x[3];
dxdt[2] = p1*x[1] - p2*x[2]*x[3] - p3*(x[2])^2;
dxdt[3] = p3*(x[2])^2;
return dxdt;
}
}
data {
int<lower=1> nsub;
int<lower=1> len[nsub];
int<lower=1> ntot;
real ts[ntot];
real obs[ntot];
}
transformed data {
int i1[nsub];
int i2[nsub];
real t0 = 0.0;
real xr[0];
int xi[0];
i1[1] = 1;
i2[1] = len[1];
for (i in 2:nsub) {
i1[i] = i2[i-1] + 1;
i2[i] = i1[i] + len[i] - 1;
}
}
parameters {
/* p1=0.04, p2=1e4, and p3=3e7 */
real<lower = 0> y0_mu;
real<lower = 0> y0_1[nsub];
real<lower = 0> sigma;
}
transformed parameters {
vector[3] y0;
vector[3] x[ntot];
vector[ntot] x3;
real theta[3] = {0.04, 1.0e4, 3.0e7};
for (i in 1:nsub) {
y0[1] = y0_1[i];
y0[2] = 0.0;
y0[3] = 0.0;
x[i1[i]:i2[i]] = ode_bdf_tol(reaction, y0, t0, ts[i1[i]:i2[i]], 1.e-8, 1.e-10, 10000, theta);
}
for (i in 1:ntot) {
x3[i] = x[i][3];
}
}
model {
y0_mu ~ lognormal(log(2.0), 0.5);
for (i in 1:nsub) {
y0_1[i] ~ lognormal(y0_mu, 0.5);
}
sigma ~ cauchy(0, 1);
obs ~ lognormal(log(x3), sigma);
}
```

data.R

```
nsub <- 8
len <- c(5, 5, 5, 5, 5, 5, 5, 5)
ntot <- 40
ts <- c(4.000e-01, 4.000e+00, 4.000e+01, 4.000e+02, 4.000e+03,
4.000e-01, 4.000e+00, 4.000e+01, 4.000e+02, 4.000e+03,
4.000e-01, 4.000e+00, 4.000e+01, 4.000e+02, 4.000e+03,
4.000e-01, 4.000e+00, 4.000e+01, 4.000e+02, 4.000e+03,
4.000e-01, 4.000e+00, 4.000e+01, 4.000e+02, 4.000e+03,
4.000e-01, 4.000e+00, 4.000e+01, 4.000e+02, 4.000e+03,
4.000e-01, 4.000e+00, 4.000e+01, 4.000e+02, 4.000e+03,
4.000e-01, 4.000e+00, 4.000e+01, 4.000e+02, 4.000e+03)
obs <- c(1.48e-02, 9.50e-02, 2.70e-01, 5.80e-01, 8.29e-01,
1.38e-02, 9.41e-02, 2.60e-01, 5.20e-01, 8.10e-01,
1.35e-02, 9.44e-02, 2.90e-01, 5.40e-01, 8.10e-01,
1.50e-02, 9.45e-02, 2.88e-01, 5.58e-01, 8.49e-01,
1.60e-02, 9.46e-02, 2.81e-01, 5.53e-01, 8.16e-01,
1.29e-02, 9.48e-02, 2.90e-01, 5.51e-01, 8.13e-01,
1.40e-02, 9.50e-02, 2.54e-01, 5.49e-01, 8.13e-01,
1.45e-02, 9.10e-02, 2.19e-01, 5.58e-01, 8.68e-01)
```

init.R

```
y0_mu <- 2.0
y0_1 <- c(2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0)
sigma <- 0.5
```

with `CVodeSetSensErrCon`

clause in stan_math

```
Elapsed Time: 322.134 seconds (Warm-up)
309.048 seconds (Sampling)
631.182 seconds (Total)
```

w/o `CVodeSetSensErrCon`

clause in stan_math

```
Elapsed Time: 225.487 seconds (Warm-up)
220.647 seconds (Sampling)
446.134 seconds (Total)
```

The performance drops ~40%.