If $f$ is complex analytic on $S=\{x+iy : |x|<1, |y|<1\}$, continuous on $\bar{S}$ and bounded by $1,2,3,4$ on each side, then is $|f(0)|>2$ possible?

I think that actually one can construct such an example using a little harmonic function theory.

First, notice that by RMT and by symmetry there is a conformal map that extends to a homeomorphism from the closed unit disc to the closed square (and which is actually conformal everywhere outside the vertices), $F:\mathbb D \to S, F(0)=0, F(\pm 1, \pm i)$ being the vertices of the square in counterclockwise order corresponding to the circle ordering (so if one fixes the image of one vertex say $F(1)=(-1,-1)$, the others are fixed eg $F(i)=(1,-1)$ etc).

(One can write a formula for this as $F(z)=c\int_0^z\frac{dz}{\sqrt{1-z^4}}$ but that is not needed)

Consider the finite non-negative measure on the unit circle given by $0, \log 2, \log 3, \log 4$ on the four open arcs $(1,i), (i-1), (-1,-i), (-i,1)$ and by zero in the four points or if you want the absolute continuous one given by $d\mu=qdt$ where $q$ takes the given values on the four open arcs and is unimportant what finite value we give it at the $4$ roots of unity of order $4$

Let $u_1(re^{i\theta})=\frac{1}{2\pi}\int_0^{2\pi}\frac{1-r^2}{1-2r\cos (\theta-t)+r^2}d\mu(t)$ the Poisson transform of $d\mu$ which is harmonic, bounded and positive in the open disc and satisfying $u_1(re^{it}) \to q(t)$ non-tangentially outside the four special points; actually note that $u_1(0)=\frac{1}{2\pi}\int_0^{2\pi}d\mu(t)=\frac {\log 2+\log 3+\log 4}{4} >2$.

So considering $g=u+iv$ holomorphic in the unit disc, $h(z)=e^{g(z)}$ almost satisfies the required properties on the four arcs since $|h(z)|=e^{u(z)}$ and $|h(0)|=e^{\frac{\log 24}{4}}=24^{1/4}>2$ (however $h$ is not continuous on the boundary, though it is "almost" so) and then clearly $f(w)=h(F^{-1}(w))=e^{g(F^{-1}(w))}$ almost satisfies the requires properties on the square and $|f(0)|>2$;

But now it is clear that $h_1(z)=(1-\epsilon)h(rz)$ for $r$ close enough to $1$ and $\epsilon>0$ small will do (will satisfy the boundness properties on the arcs and will be continuous, even holomorphic on the closed disc) and then $|h_1(0)|>2$ if $1-r$ hence $\epsilon$ are small enough, so taking $f_1(w)=h_1(F^{-1}(w))$ solves the problem.

Note that since $F, F^{-1}$ are not conformal at the vertices, $f_1$ is only continuous on the closed square (though holomorphic outside the vertices) despite that $h_1$ is holomorphic on the closed unit disc