The ODE $y'(x)+P(x)y(x)=Q(x)$ has solution $$I(x)y(x)=\int I(x)Q(x)\,dx$$ where $I(x)=\exp\int P(x)\,dx$. Equivalently, $$Y(x)+P(x)\int_0^xY(t)\,dt=Q(x)\tag1$$ has solution $$Y(x)=\frac d{dx}\frac{\int I(x)Q^*(x)\,dx}{I(x)}=\frac{I(x)^2Q^*(x)-I'(x)\int I(x)Q^*(x)\,dx}{I(x)^2}$$ where $Y=y'$ and $Q^*(x)=Q(x)+P(x)y(0)$. Equation $(1)$ gives the limiting case, where $$\int_0^xY(t)\,dt=\lim_{n\to\infty}\frac xn\sum_{k=0}^nY\left(\frac{kx}n\right).$$ Given $P(x),Q(x)$, what could be said about the solutions of the functional equation $$Y(nx)+\frac{xP(x)}n\sum_{k=0}^nY(kx)=Q(x)\tag2,$$ where $n$ is no longer under the limit?
That is, what is the behaviour of the families of solutions to $(2)$ as $n$ increases?
(Cross-posted on MathOverflow.)