I would like to compute the solution of an ODE at several (irregularly spaced) times `ts`

. This computation can currently be done in a loop:

```
for (t in ts)
y[t] = matrix_exp(t * A) * y0;
```

My question is: Would it be possible to reduce the computational load by vectorizing this computation? I assume @yizhang is probably the most knowledgeable on this topic.

This could naturally be included in the `scale_matrix_exp_multiply`

function that recently was exposed to Stan, i.e. it would get the signature:

```
scale_matrix_exp_multiply(vector t, matrix A, matrix B)
```

If it is expected that this can reduce computationally load, I can submit a feature request for this.