The solutions of $(z-i)^n+(z+i)^n=z^n$ are real.

Define: $$ g(x) := \cos\left(n\cdot\arccos\frac{x}{\sqrt{x^2+1}}\right) $$ We claim that we can partition $\mathbb{R}$ into $n$ intervals, where in each interval $g$ continuously moves between the values $-1$ and $1$. To show this, let $f(x):=\frac{x}{\sqrt{x^2+1}}$. Note that $f$ increases from $-1$ to $1$ as $x$ goes from $-\infty$ to $\infty$. This means that $n\cdot\arccos\left(f(x)\right)$ decreases from $n\pi$ to $0$, as as $x$ goes from $-\infty$ to $\infty$, which proves our claim.

Next, define: $$ h(x) := \frac{1}{2}\left(\frac{x}{\sqrt{x^2+1}}\right)^{n} $$ We note that $-\frac{1}{2}\lt h\lt \frac{1}{2}$ for all real $x$. We use the intermediate value theorem to conclude that $h$ intersects $g$ in each of the $n$ intervals metioned in the previous paragraph, i.e. $g(x)=h(x)$ has $n$ real solutions.

Let $x_0\in\mathbb{R}$ be such that $h(x_0)=g(x_0)$. Then, by setting $\theta :=\arccos\frac{x_0}{\sqrt{x_0^2+1}}$, we get: $$ \begin{equation} \begin{split} (x_0+i)^n+(x_0-i)^n &=(x_0+i)^n+\overline{(x_0+i)^n} \\ \\ &=2\cdot\Re\left((x_0+i)^n\right) \\ \\&=2\cdot\Re\left(e^{in\theta}\left(\sqrt{x_0^2+1}\right)^n\right)\\ \\ &= 2\cos\left(n\theta\right)\left(\sqrt{x_0^2+1}\right)^n\\ \\ &=x_0^n \end{split} \end{equation} $$ This implies that the original equation has at least $n$ real solutions. However, since it is a polynomial of degree $n$, it cannot have more than $n$ solutions, which means that all solutions must be real.

Consider the function $p_n:\Bbb{R}\to\Bbb{C}$ given by $$ p_n(x)=(1+ix)^n. $$ The argument $$ \theta_n(x):=\arg p_n(x)=n\arctan x $$ then increases from $-n\pi/2$ to $+n\pi/2$ as $x$ increases from $-\infty$ to $+\infty$.

Observation #1: In any interval, where $\theta_n(x)$ either

• increases from $-\pi/2+2k\pi$ to $-\pi/3+2k\pi$ for some integer $k$, or
• increases from $2k\pi+\pi/3$ to $2k\pi+\pi/2$,

there exists an $x$ such that the real part of $p_n(x)=1/2$.

Proof: The absolute value $p_n(x)$ is always $\ge1$. So if $[a,b]$ is an interval of one of the above types, then at the end corresponding to an odd multiple of $\pi/2$, the real part of $p_n(x)=0$. At the other end the cosine of the argument is $1/2$ and therefore the real part of $p_n(x)\ge1/2$. The claim follows from continuity.

Observation #2: If $p_n(x)$ has real part $=1/2$, then $z=1/x$ is a solution of the original equation.

Proof: If $p_n(x)=1/2+it$, then $$ (1+ix)^n+(1-ix)^n=p_n(x)+\overline{p_n}(x)=(\frac12+it)+(\frac12-it)=1. $$ Dividing this by $x^n$ gives $$ (1/x+i)^n+(1/x-i)^n=(1/x)^n $$ proving the claim.

The claim follows from these observations, the known codomain of $\theta_n(x)=n\arctan x$, and the fact that $z=0$ is a solution for all odd $n$. If $n=2m$ is even, then there are exactly $2m$ intervals of the type covered in the first observation. OTOH, if $n=2m+1$ then those same $2m$ intervals are included in the range of $\theta_n(x)$, and we have the extra zero from $z=0$.

In either case we have found $n$ real zeros, so there cannot be others.