From old qualifying exam: Let $E$ be the union of the two coordinate axes, i.e. $E = \{z=x+iy : xy=0\}$. Describe all entire functions satisfying $f(E) \subset E$.
I feel like the best approach is to consider the power series of $f$. My first approach was to write down constraints by considering the function applied to the real or imaginary axis. When I didn't get anywhere with this, I began thinking of the function geometrically: on $E$ it's allowed to scale by a real constant, and rotate by $k\pi/2$. But again, I couldn't see how to usefully translate this to produce information about the power series. Thanks!
As an example, $z^2$ has this property. In fact, so does $az^2+bz^4$ (with $a,b \in \mathbb{R}$) since each term maps the imaginary axis to the real axis, which ends up back on the real axis when added. A similar argument shows real odd polynomials work, too.
If $f(E)\subset E$, then $(f^2)(\mathbb{R})\subset\mathbb{R}$. Thus, the coefficients of the power series of $(f)^2$ (around $z=0$) are real. If $\{c_n\}$ are the Taylor coefficients of $f$ around $z=0$, then $$ \sum_{j=0}^nc_jc_{n-j}\in\mathbb{R},\quad n=0,1,2,\dots\tag1 $$ If $f(E)\subset E$, the same is true of $g(z)=f(i\,z)$. Since the coefficients of $g$ are $i^n\,c_n$, it folows that $$ i^n\,\sum_{j=0}^nc_jc_{n-j}\in\mathbb{R},\quad n=0,1,2,\dots\tag2 $$ From (1) and (2) we get $$ \sum_{j=0}^nc_jc_{n-j}=0\quad\text{if $n$ is odd.}\tag3 $$ Now it is not to difficult to see that the power series of $f$ must be one of $$ \sum_{n=0}^\infty c_nz^{2n},\quad i\,\sum_{n=0}^\infty c_nz^{2n},\quad \sum_{n=0}^\infty c_nz^{2n+1},\quad i\,\sum_{n=0}^\infty c_nz^{2n+1},\text{ with } c_n\in\mathbb{R}. $$