I have an ODE system defining a mathematical model of a biological system, say
$$ \frac{da}{dt}=f_1(a,b,\ldots,z,p)\\ \frac{db}{dt}=f_2(a,b,\ldots,z,p)\\ \cdots\\ \frac{dz}{dt}=f_n(a,b,\ldots,z,p) $$
with state variables, $a,b,\ldots,z$, and parameter vector, $p$.
In the end, I need to calculate a scalar model response, $f$, defined as
$$ f(p)=\int_0^{t_\text{end}}dt \left(a(t)+b(t)\right)+\int_0^{t_\text{end}}dt \left(x(t)+z(t)\right) $$
Question: what is the best way to calculate $f(p)$ apart from numerical integration once the time-courses are calculated?