$\newcommand\de\delta$Let us also present the Fourier transform argument, assuming that $|F|$ is bounded by a polynomial, so that $F$ may be considered a (tempered) distribution (in the generalized-function sense). Let then $\hat F$ denote the Fourier transform of $F$.
Equation (1) in the other answer yields $$c_n\de(t)=\hat G_n(t)=e^{ita_n}\hat F(t)+e^{-ita_n}\hat F(t)-2\hat F(t) =2\hat F(t)(\cos ta_n-1),$$ where $\de$ is the delta function.
If the equality $\cos ta_n-1=0$ takes place for some real $t$ and all $n$, then $t=0$ (since the $a_n$'s are nonzero and go to $0$). So, the support of $\hat F$ is $\{0\}$. So (see e.g. "For every compact subset $K\subseteq U$ there exist constants $C_{K}>0$ and $N_{K}\in \mathbb {N}$ such that for all $f\in C_{c}^{\infty }(U)$ with support contained in $K$ [...]" here), we have $\hat F=\sum_{j=0}^n a_j\de^{(j)}$ for some $n\in\{0,1,\dots\}$ and some complex $a_j$'s, where $\de^{(j)}$ is the $j$th derivative of the delta function $\de$. So, $F$ is a polynomial. Since the second difference $G_n$ of $F$ is (the) constant ($c_n$), it follows that the polynomial $F$ is quadratic. Thus, $f=F'$ is an affine function. (Vice versa, any affine function $f$ satisfies your system of functional equations.)