This is the first draft of allowing `ts`

in ODE solvers to be parameters. The motivation is to make inference on the time of observations. What I did is to attach the RHS of the ODE to the `decouple_states`

methods as that’s the sensitivity of the ODE variables with respect to the time.

The test model below shows what the implementation tires to achieve

```
functions {
real[] sho(real t,
real[] y,
real[] theta,
real[] x,
int[] x_int) {
real dydt[2];
dydt[1] = y[2];
dydt[2] = -y[1] - theta[1] * y[2];
return dydt;
}
}
data {
int<lower=1> T;
real y0_d[2];
real t0;
real theta_d[1];
real y_hat[T];
real x[0];
int x_int[0];
}
parameters {
real<lower=0, upper=0.2> ts1;
real<lower=ts1, upper=0.3> ts2;
real<lower=ts2, upper=0.5> ts3;
}
transformed parameters{
real<lower = 0> ts[T];
for(i in 1:T) ts[i] = 0.1*i;
ts[1] = ts1;
ts[2] = ts2;
ts[3] = ts3;
}
model {
real sol[T,2];
sol = integrate_ode_rk45(sho, y0_d, t0, ts, theta_d, x, x_int);
y_hat ~ normal(sol[1:T, 2], 0.01);
}
```

With the ODE solution data, parameters `ts1`

, `ts2`

, and `ts3`

are inference results of the first three time steps.

Comments and suggestions?