Consider the following first order linear ODE: $$\vec{x}'(t) = A(t)\vec{x}(t)+\vec{f}(t).$$ Let $\Psi(t)$ be a fundamental matrix, i.e. a matrix whose columns are linearly independent solutions to the above equation when $\vec{f} \equiv \vec{0}$.
By making the ansazt $\vec{x}(t)=\Psi(t)\vec{c}(t)$ for the solution subject to the initial condition $\vec{x}(t_0)=\vec{x}_0$, it is trivial to show that we get $$\vec{x}(t)=\Xi(t,t_0)\vec{x}_0+\int_{t_0}^{t}\Xi(t,\tau)\vec{f}(\tau)d\tau$$ where we have defined $\Xi(t,\tau)=\Psi(t)\Psi^{-1}(\tau)$.
I now want to consider the linear recurrence analog $$\vec{x}(k+1) = A(k)\vec{x}(k)+\vec{f}(k). $$ Again, let $\Psi(k)$ be a fundamental matrix, i.e. a matrix whose columns are linearly independent solutions to the above equation when $\vec{f} \equiv \vec{0}$.
It appears natural to make the ansatz $\vec{x}(k) =\Psi(k)\vec{c}(k)$. We note that $\Psi(k+1)=A(k)\Psi(k)$, so we get $\vec{x}(k+1) = \Psi(k+1)\vec{c}(k+1) = A(k)\vec{x}(k)+A(k)\Psi(k)[\vec{c}(k+1)-\vec{c}(k)] = A(k)\vec{x}(k)+\vec{f}(k)$, whence $\vec{c}(k+1) = \vec{c}(k)+\Psi(k)^{-1}A(k)^{-1}\vec{f}(k).$
However, this doesn't seem quite right since we shouldn't have to assume that $A(k)$ is invertible.
My guess is that the final formula will look something like $$\vec{x}(n) = \Xi(n,n_0)\vec{x}_0+\sum_{\eta = n_0}^{n-1}\Xi(n,\eta)\vec{f}(\tau)d\tau$$ if we prescribe the initial condition $\vec{x}(n_0)=\vec{x}_0$. Here we should have $\Xi(n+1,\eta) = A(n)\Xi(n,\eta)$ as well as $\Xi(n,n)=I$. Indeed, with these properties, we get $\vec{x}(n+1)=A(n)\vec{x}(n)+\vec{f}(n)$. In particular, if $A(n)=A$ is constant, we have $\Xi(n,\eta)=A^{n-\eta}, \ n \geq \eta$.
However, I am not sure hot to derive it. Any suggestions or comments?
Notice that in the smooth case the expression satisfies the product rule:
$$x(t) = \Psi(t)c(t)\implies \dot x(t) = \dot\Psi(t) c(t) + \Psi(t)\dot c(t)$$
However with the ansatz $x(k) = \Psi(k)c(k)$ we simply have $x(k+1) = \Psi(k+1)c(k+1)$ from which the problems you observed stem from. I one ougth to work with a finite difference operator here instead. But let's just see what happens if we write down the first few iterations
Which suggests that $$x_k = \Big[\prod_{i=k-1}^{0}A_i\Big] x_0 + \sum_{i=0}^{k-1}\Big[\prod_{j=k-1}^{i+1} A_j\Big] f_i$$
which is in line with your guess. Indeed if we introduce $\Xi(k,l) = A_{k-1} \cdot \ldots \cdot A_{l}$ for $l\le k-1$ we can write the solution as
$$x_k = \Xi(k,0)x_0 + \sum_{i=0}^{k-1} \Xi(k,i+1) f_i $$
Note that indeed we have $\Xi(k+1,l) = A_k \Xi(k,l)$ and $\Xi(k,k) = I$ if we interpret this case as an empty product. Thereupon the proof is a very simple induction: \begin{align} x_{k+1} = A_k x_k + f_k &= A_k \Xi(k,0)x_0 + f_k + \sum_{i=0}^{k-1} A_k\Xi(k,i+1) f_i \\ &= \Xi(k+1,0)x_0 + f_k + \sum_{i=0}^{k-1}\Xi(k+1,i+1) f_i \\ &= \Xi(k+1,0) + \sum_{i=0}^{k}\Xi(k+1,i+1) f_i \\ \end{align}
You can find more details in this script from Stanford University