[update - the alternative strategy I propose in the concluding remark seems not to work, as I explain below, though the main proof itself is fine.]
As you have observed, taking derivatives of each side (via the second fundamental theorem per Sarvesh Ravichandran Iyer's comment) gives us $f'(x)=f(x+1)-f(x)$. Note this implies $f'$ is continuous, and since $|f(x)|\leq M$ for some $M$ we also have $|f'(x)|\leq 2M$. In particular, $f$ is $2M$-Lipschitz.
One has to argue a little more carefully to show the Taylor series actually converges to $f(x)$ for each $x$. This can be done, as previously noted by Gerd, via Taylor remainder estimates, using the bound $|f^{(k)}(x)|\leq 2^kM$. However, as mentioned in the comments, this only immediately shows you that $\sum_{k=2}^\infty \frac{f^{(k)}(x)}{k!}=0$. It is not clear why each term must be $0$. Therefore, while your approach may still work, it is not obvious to me how to salvage it.
However, just the bound on the first derivative is already enough to prove the result, without taking further derivatives. Here is a proof:
Proof
Our strategy will be to show that if $f'$ is ever (WLOG) positive, then we can find another point where $f'$ is quantitatively larger, so we obtain a sequence of points on which $f'$ increases without bound, contradicting the bound $|f'(x)|\leq 2M$.
Suppose by contradiction that WLOG $f'(x)=\epsilon>0$ for some $x$. Then for $t\in[x+1-\frac{\epsilon}{2M},x+1]$, we have by the Lipschitz condition that
$$f(t)\geq f(x+1)-2M(x+1-t),$$
which is linear from $f(x)$ to $f(x+1)$ on that interval,
and so we obtain $$\int_{x+1-\frac{\epsilon}{2M}}^{x+1} (f(t)-f(x))dt \geq \frac{1}{2}\cdot \frac{\epsilon}{2M}(f(x+1)-f(x))=\frac{\epsilon^2}{4M}.$$
Since the average of $f$ over $[x,x+1]$ is $f(x)$, we must have
$$\int_x^{x+1-\frac{\epsilon}{2M}} (f(x)-f(t)) dt\geq \frac{\epsilon^2}{4M}.\tag{1}$$
Now, if $y$ is the largest value in $[x,x+1]$ for which $f(y)=f(x)$, then we obtain $$\int_x^y (f(x)-f(t)) dt\geq \frac{\epsilon^2}{4M},\tag{2}$$
since the integrand in (1) is negative after $t=y$.
However, the Lipschitz condition gives us $$f(x)-f(t)\leq 2M(t-x),$$ so integrating the right hand side over $[x,y]$ and combining with (2) we have the estimate
$$M(y-x)^2\geq \frac{\epsilon^2}{4M},$$ so that $y-x\geq \frac{\epsilon}{2M}$.
Then by the mean value theorem we have some $z\in [y,x+1]$ with
$$f'(z)=\frac{f(x+1)-f(y)}{x+1-y}\geq\frac{\epsilon}{1-\frac{\epsilon}{2M}}.$$
Note also that if we started with $f'(x)=\epsilon'\geq\epsilon$, then we would conclude $$f'(z)\geq \frac{\epsilon'}{1-\frac{\epsilon'}{2M}}\geq \frac{\epsilon}{1-\frac{\epsilon}{2M}},$$ so we can summarize that if there is a point $x$ with $f'(x)\geq \epsilon$, then there is a point $z$ with $f'(z)\geq \frac{\epsilon}{1-\frac{\epsilon}{2M}}$.
Then letting $\epsilon_0=\epsilon$ and $\epsilon_{n+1}=\frac{\epsilon_n}{1-\frac{\epsilon_n}{2M}}$, we have for each $n$ a point $x_n$ with $f'(x_n)\geq \epsilon_n$. But $\epsilon_n$ increases without bound - to see this note that $\epsilon_n$ is certainly increasing, and so we may estimate
$$\frac{\epsilon_{n+1}}{\epsilon_n}= \frac{1}{1-\frac{\epsilon_n}{2M}}\geq \frac{1}{1-\frac{\epsilon}{2M}}>1,$$
so $$f'(x_n)\geq \epsilon_n\geq \epsilon\left(\frac{1}{1-\frac{\epsilon}{2M}}\right)^n\to\infty.$$
This contradicts the boundedness of $f'(x)$.
Remark
Another avenue of proof might be to observe that if $\mathcal F\subseteq C_b(\mathbb R)$ is the space of bounded continuous functions satisfying the given condition, then $\mathcal F$ is closed under the supremum norm, and the differential $f\mapsto f'$ restricts to a continuous (in the sup norm) linear map $D\colon \mathcal F\to \mathcal F$, so by this result (*not quite actually), $\mathcal F$ is finite dimensional.
[*edit: just realized that at the link they are discussing functions in $C([0,1])$, so the result might not generalize to $C_b(\mathbb R)$, I’ll have to think more on that. In fact, I think this is false, and you should actually be able to embed an isometric copy of $c_b$ (space of bounded sequences) into $C_b(\mathbb R)$ so that the image consists of smooth functions and differentiation is continuous. Thus this alternative strategy is probably not salvageable, but I’ll leave the idea here anyway for just a bit before deleting, in case I or anyone else thinks of a fix. We might be able to get some mileage out of the fact that $\mathcal F$ consists of real analytic functions, and hence the restriction to a compact interval is injective.]
Then if $P$ is the characteristic polynomial of $D$, we have $P(D)=0$, so that every member of $\mathcal F$ lies in the solution set to a homogenous linear $n$-th order ODE with constant coefficients. Since the functions in $\mathcal F$ are bounded, we then know that $\mathcal F$ consists entirely of functions of the form $$f(x)=C+\sum_i a_i\sin(\omega_ix)+b_i\cos(\omega_i x).$$
One then simply needs to show that any such functions satisfying the given condition are constant. Ex post facto, this is true, as we have proven the result with other means, but I didn't immediately see how to cleanly show this (I doubt it is hard, but it might or might not be messy) and may not have time to further explore, so I will leave this as a remark.
