Such a function does not exist. Assume first that $f$ is real on the real line. We use the theorem of Marchenko–Ostrovski that if $f$ is entire and all zeros and solutions of $f(z)=3$ are real, then $f$ belongs to the Laguerre–Pólya class (closure of polynomials with all zeros real).
This theorem also tells us how the graph of this function on the real line looks like: if
$a_k$ are zeros and $b_k$ are 3-points, then
their order on the real line is
$$\dotsb a_k\leq a_{k+1}<b_k\leq b_{k+1}< a_{k+2}\leq\dotsb$$
with some choice of enumeration, and there is one critical point in each interval $[a_k,a_{k+1}]$ and in $[b_k,b_{k+1}]$, and no other critical points. The sequences of zeros and $3$-points are infinite in both directions.
In particular,
if $f(n)=0$, and $f'(n)\geq 0$ then $f(n+k)=k, 0\leq k\leq 3$, $f(n+4)=3$, $f(n+4+k)=3-k,\; 0\leq k\leq 3$ and so on.
Now fix $n$ such that $f(n)=0$, and consider the function on the interval $[n,n+3]$. Since $f(n+k)=0$, $0\leq k\leq 3$, the function $g(x)=f(x)-x+n$ has $4$ zeros on the interval, therefore $g'$ has at least $3$ zeros, and $g''=f''$ has at least $2$ zeros. This contradicts the fact that for functions of Laguarre–Pólya class, the zeros of $f^{(n)}$ and $f^{(n+1)}$ are interlacent: the adjacent zeros of $f'$ lie in $[n-1,n]$ and $[n+3,n+4]$,
and thus $f''$ can have at most one zero on $[n,n+3]$.
Reference for the Marchenko–Ostrovski theorem:
Zbl 327.34021
Marchenko, V. A.; Ostrovskii, I. V.
A characterization of the spectrum of Hill’s operator. (Russian)
Mat. Sb., N. Ser. 97(139), 540-606 (1975).
(English translation).
Remark 1. Marchenko–Ostrovski theorem has an a priori assumption that $f$ maps the real line into itself. However this property follows from our assumption that $f$-preimage of $4$ points belongs to the real line, see, for example,
Zbl 1172.30004
Bergweiler, Walter, Eremenko, Alexandre
Meromorphic functions with linearly distributed values and Julia sets of rational functions.
Proc. Amer. Math. Soc. 137 (2009), no. 7, 2329–2333. arXiv.
where it is proved that whenever all solutions of
$f(z)=a$ are real, for $4$ distinct values of $a$,
the function maps the real line to itself.
Remark 3. Analysis of the above proof shows that a slightly stronger result holds: If $f(\mathbf{Z})\subset\mathbf{Z}$ and $f^{-1}(\{0,1,2,3\})\subset\mathbf{Z}$ then $f(z)=z+b$,
where $b\in\mathbf{Z}$.