I have the following differential equation in the form $\mathbf{y}'=\mathbf{A}(t)\mathbf{y}$:
$$ \frac{d}{dt} \begin{bmatrix} y_0 \\ y_1 \\ y_2 \\ \vdots \\ y_n \end{bmatrix} = \begin{bmatrix} -\lambda_s -\lambda_b[f_1(t) + \cdots + f_n(t)] & 0 & 0 & \cdots & 0 \\ \lambda_bf_1(t) & 0 & 0 & \cdots & 0 \\ \lambda_bf_2(t) & 0 & 0 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & 0\\ \lambda_bf_n(t) & 0 & 0 & 0 & 0 \end{bmatrix} \begin{bmatrix} y_0 \\ y_1 \\ y_2 \\ \vdots \\ y_n \end{bmatrix} $$
Where $\lambda_s$ and $\lambda_b$ are constants and $f_i(t) \in [0, 1]$ are functions defined for all $t$.
I want to calculate $\mathbf{y}(t=T)$, where $T > 0$, solely with the following information:
- The constants $\lambda_s$, $\lambda_b$, $T$, and the initial condition $\mathbf{y}(t=0)$
- The integrals $\int_0^T f_i(t)\ dt$ or any other integral in the form $\int_0^T g_i(\mathbf{y}, t)\ dt$ but critically, $\lambda_b$ and $\lambda_s$ CANNOT be in the expression. This is because the purpose of all of this is to see how $\mathbf{y}(t=T)$ varies as a function of $\lambda_b$ and $\lambda_s$.
For $y_0$, the differential equation is separable and it is easy to find $y_0(T)$ with only the above information:
$$ \frac{dy_0}{dt} = -y_0\left(\lambda _s + \lambda_b \sum_{i=1}^{n} f_i(t)\right) $$
$$ y_0(T) = y_0(t=0)\cdot \exp{\left(-\lambda_s T + \lambda_b\sum_{i=1}^n \int_0^T f_i(t)\ dt \right)} $$
If there wasn't the requirement for $\lambda_b$ and $\lambda_s$ to not be present in the integrals, then the following solution would exist for $k \in [1, n]$:
$$ y_k(T) = y_0(t=0)\cdot \lambda_b \int_0^T f_k(t)\cdot \exp{\left(-\lambda_s T + \lambda_b\sum_{i=1}^n \int_0^t f_i(\tau)\ d\tau \right)}\ dt + y_k(t=0) $$
But as it stands, this solution does not help because I can't evaluate this for varying $\lambda$ with the integrals pre-computed. The fact that this form of solution exists for $y_0$ makes me think it might be possible for the other terms, and according to this Wikipedia article, a general solution
$$ \mathbf{y}(T)=\exp{\left(\int_0^T \mathbf{A}(t)\ dt\right)}\mathbf{y}(t=0) $$
exists, but the solution to this is different to the numerical solution of a test case I tried. I wonder if it's because $\mathbf{A}$ is singular, but I really have no idea.
The reason for these constraints is that I want to calculate $\mathbf{y}(T)$ for billions of vectors $\mathbf{y}$ and therefore can only store the minimum number of variables for each $\mathbf{y}$ (e.g. integral quantities).