Is there a simple way to prove this fact :
For all holomorphic functions $f : \mathbb C \to \mathbb C$, there is an holomorphic function $\psi : \mathbb C \to \mathbb C$ such that $$\psi(z+1) = \psi(z) + f(z) $$
The solution I know use Galois covering spaces and summand of automorphy.
Thanks in advance.
As I mentioned in a comment, there are two very simple ways of solving the functional equation $\psi(z+1)=\psi(z)+f(z)$, for $\psi$ in terms of $f$. They are, \begin{align} \psi(z)&=-f(z)-f(z+1)-f(z+2)-f(z+3)-\cdots,{\rm\ and,}\\ \psi(z)&=f(z-1)+f(z-2)+f(z-3)+\cdots. \end{align} These solutions will only apply if either of the respective sums converge and, if we want $\psi$ to be holomorphic, then we require locally uniform convergence. For example, taking $f(z)=e^{az}$ then, for $\Re(a)\not=0$, the solutions above give $\psi(z)=e^{az}/(e^a-1)$ (using the first method for $\Re(a)<0$ and the second for $\Re(a)>0$).
More generally, we can prove the following,
Lemma: Let $f\colon\mathbb{C}\to\mathbb{C}$ be entire. Then, there exists entire $g,h$ with $f=g+h$ and such that, for each $K>0$, $z^2g(z)$ and $z^2h(-z)$ are uniformly bounded on the horizontal strip $-K\le\Im(z)\le K$ with $\Re(z)>0$.
I'll prove this in a moment. First, however, we can see that it gives an easy solution to the functional equation by applying the first method above to $g$ and the second to $h$, \begin{align} \psi(z)=&-g(z)-g(z+1)-g(z+2)-\cdots\\ &+h(z-1)+h(z-2)+h(z-3)+\cdots. \end{align}
I'll now prove the lemma stated above, using a result on approximating functions by entire functions due to Arakelian. Let $E$ be the set of $z\in\mathbb{C}$ such that $2\lvert\Im(z)\rvert+1\le\lvert\Re(z)\rvert$. This is an unbounded closed subset of $\mathbb{C}$ and consists of two connected components, one for $\Re(z)\le-1$ and one for $\Re(z)\ge1$. Furthermore, for each $K>0$, the strip $-K\le\Im(z)\le K$ is contained in the union of $E$ and a bounded set. So, it is enough that the functions $g,h$ constructed in the lemma are such that $z^2g(z)$ and $z^2h(-z)$ are bounded on the component of $E$ with $\Re(z)\ge1$.
We can define a continuous function $F\colon E\to\mathbb{C}$ by setting $F(z)=e^{z^2}f(z)$ for $\Re(z)\le-1$ and $F(z)=0$ otherwise. This is also analytic on the interior of $E$.
Now, $E$ is an Arakelian set. That is, it is a closed subset of $\mathbb{C}$ whose complement does not contain any bounded connected components and such that, for each closed disk $D$, the union of the closed bounded components of $\mathbb{C}\setminus(E\cup D)$ is bounded (actually, it is empty, so trivially bounded). Then, Arakelian's approximation theorem states that, for any $\epsilon > 0$, there is an entire function $G$ with $\lvert F-G\rvert\le\epsilon$ on $E$. If we set $g(z)=e^{-z^2}G(z)$ and $h(z)=f(z)-g(z)$ then we see that $\lvert e^{z^2}g(z)\rvert$ and $\rvert e^{z^2} h(-z)\rvert$ are bounded by $\epsilon$ on the component of $E$ with $\Re(z)\ge1$.
Finally, on $E$, $$ \Re(z^2)=\Re(z)^2-\Im(z)^2\ge\frac35\left(\Re(z)^2+\Im(z)^2\right)=\frac35\lvert z\rvert^2. $$ So, $$ \lvert e^{z^2}\rvert=\exp(\Re(z^2))\ge\exp(3\lvert z\rvert^2/5)\ge\frac35\lvert z\rvert^2, $$ and $z^2g(z)$, $z^2h(-z)$ are bounded on the component of $E$ with $\Re(z)\ge1$ as required.
(This is nothing but an expansion of my comment; it clarifies what can easily be done by a naïve power series method)
The question is to prove that the finite difference operator $\Delta : f \mapsto f(z+1) - f(z)$ defines a surjection $\mathcal O(\mathbb C) \to \mathcal O(\mathbb C)$. ($\mathcal O(\mathbb C)$ is the algebra of entire functions, that is holomorphic functions on all of $\mathbb C$).
Because of standard complex analysis theorems, an entire function is the sum of a power series of infinite convergence radius. That is, series summation gives an isomorphism $$ A_\infty = \left\{\sum_{n\geq 0} a_n z^n \in \mathbb C[[z]] \, \middle | \, \forall R > 0, a_n = O(R^{-n}) \right \} \to \mathcal O(\mathbb C).$$
So the finite difference operator also defines an operator $\Delta : A_\infty \to A_\infty$. Basically, it is given by Newton's binomial formula: for $j \geq 0$, we have $$\Delta(z^j) = (z+1)^j - z^j = \sum_{i=0}^{j-1} \binom j i z^i.$$ For a power series $\sum_{j \geq 0} a_j z^j \in \mathbb C[[z]]$, one can apply the last formula and linearity to get $$\Delta\left( \sum_{j=0}^K a_j z^j \right) = \sum_{j=0}^K a_j \left( \sum_{i=0}^{j-1} \binom j i z^i \right) = \sum_{i=0}^{K-1} \left(\sum_{j=i+1}^K a_j \binom j i \right) z^i.$$ So, as, $K$ grows, all of the coefficients of this polynomial keep changing. Because $\binom j i$ is a degree-$i$ polynomial in $j$, the series $\sum_{j>i} a_j \binom j i$ converges and the operator $\Delta : A_\infty \to A_\infty$ is defined by the formula: $$\begin{array}{ccc} A_\infty&\to& A_\infty\\ \sum_{n \geq 0} a_n z^n & \mapsto &\displaystyle \sum_{n\geq 0} \left( \sum_{k=n+1}^{+\infty} a_k \binom k n \right) z^n.\end{array}$$
So the original question is really equivalent to:
Question: is this operator $\Delta : A_\infty \to A_\infty$ surjective?
To address this question, it seems quite reasonable to look at polynomials of degree $N$ first, then to look at $N \to \infty$. Exactly like in the computation of $\Delta$, you need to have some kind of convergence to extend the definitions from the polynomial case to the power series case. In the previous computation, we needed the convergence of $\sum_{k \geq n} a_k \binom k n$ for all $n$ and the fact that $\sum_{n \geq 0} a_n z^n$ has an infinite convergence radius was more than enough for that.
What I'll explain is that in order to get a map $I : A_\infty \to A_\infty$ such that $\Delta \circ I = \mathrm{id}_{A_\infty}$ (which is equivalent to proving that $\Delta$ is surjective) a similar approach works, but only if we can assume that the decay of the $(a_n)$ is stronger than the $\forall R > 0, a_n = O(R^{-n})$ hypothesis. That's why this isn't really an answer to the original question.
Obviously (either because of the definition or because of the big messy formula above) $\Delta$ induces a linear map $\mathbb C_{n}[z] \to \mathbb C_{n-1}[z]$. It's quite easy to show that the only polynomials in $\ker \Delta$ are the constant ones, so a dimension computation shows that $\Delta : \mathbb C_{n}[z] \to \mathbb C_{n-1}[z]$ is onto, for all $n$. There is then a map $I : \mathbb C_{n-1}[z] \to \mathbb C_{n}[z]$ such that $\Delta \circ I = \mathrm{id}_{\mathbb C_{n-1}[z]}$ and this map is unique, up to adding constant polynomials.
In fact, one can be quite explicit about that: if $B_n$ is the $n$-th Bernoulli number, the polynomial $$P_k = \frac{1}{k+1} \left( \sum_{n=0}^{k+1} \binom{k+1}{n} B_{k+1-n} z^n\right)$$ (which, up to the $1/(k+1)$ factor, is the $(k+1)$-th Bernoulli polynomial) satisfies the following relation $$\Delta P_k = z^k,$$ so it seems like a good idea to define $I(z^k) = P_k$.
If $\sum_{k\geq 0} a_k z^k \in A_\infty$, we then have that $$\begin{align*}I\left(\sum_{k=0}^K a_k z^k\right) &= \sum_{k=0}^K a_k \left(\frac{1}{k+1} \sum_{n=0}^{k+1} \binom{k+1}{n} B_{k+1-n} z^n\right)\\ & = \sum_{n=0}^{K+1} \left(\sum_{k=n-1}^K \frac{a_k}{k+1} \binom{k+1}n B_{k+1-n} \right) z^n.\end{align*}$$
Now, if the decay of $(a_k)$ was strong enough to guarantee the absolute convergence of
$\sum_{k\geq n-1} \frac{a_k}{k+1} \binom{k+1}n B_{k+1-n}$, the formula above (with $K$ replaced by $\infty$) would define a linear map $I : A_\infty \to A_\infty$ such that $\Delta \circ I = \mathrm{id}_{A_\infty}$. This would be the case if we had $|B_p| \leq C\cdot R^p$ for some constants $C$, $R$, for instance. The problem is that the real asymptotics for the Bernoulli numbers is
$$|B_{2n}| \sim 4\sqrt{\pi n}\left(\frac{n}{\pi e}\right)^{2n}$$ (for $n > 0$, $B_{2n+1} = 0$).
It is then well possible to find sequences $(a_n)$ such that $a_n = O(R^{-n})$ for all $R$, but which nevertheless satisfy $$\limsup_{k\to\infty} \left|\frac{a_k}{k+1} \binom{k+1}n B_{k+1-n}\right| = +\infty.$$ For example, $a_k = k^{-\alpha k}$ seems to work for $0 < \alpha < 1$.
So I think this method works if we restrict ourselves to classes of entire functions with strongly decaying Taylor coefficients (strongly enough for our $I$ to be defined), but not for any entire function.
Like I said in my comment, I'm surprised that the answer to the Question asked above could be "yes" with this natural method failing. If the OP hadn't claimed the original result to be true, I would by now think the $\Delta$ operator not to be surjective on $A_\infty$ or $\mathcal O(\mathbb C)$.
On the other hand, if you try this method (using polynomials of higher and higher degree to approximate your power series solution) to construct elements of $\ker \Delta$, you would probably only find the constant functions, whereas the kernel of $\Delta$ in $\mathcal O(\mathbb C)$ [or $A_\infty$] is much larger than that (for instance it cointains $z \mapsto e^{2i\pi z}$). So perhaps there is a good reason for which this method was doomed to fail (some instability property?) and my surprise is really nothing but a consequence of my naïveté.
