# calculating integral for an ODE system

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?

The common way of tracking such properties is indeed to add an additional ODE to the system, here: $$\frac{df}{dt}=a(t)+b(t)+x(t)+z(t)$$ with initial condition $$f(0)=0.$$ If this is not possible, then a smart choice of output times e.g. quadrature points in t for an integration afterwards is recommend.