How to solve $e^x=x$?

$$f(z)= z\cdot e^{-z}$$ is an entire function, hence the great Picard's theorem gives that $f(z)$ takes any complex value (with at most one exception) infinitely often. The exceptional value is obviously zero, since $f(z)=0$ implies $z=0$, so there is an infinite number of complex numbers for which $$ f(z) = 1,$$ i.e. $$ z = e^z. $$ Cauchy's integral formula gives also that the number of zeroes of $z-e^z$ inside a region bounded by a simple closed curve $\gamma$, counter-clockwise oriented, is given by: $$ N_{\gamma}=\frac{1}{2\pi i}\oint_{\gamma}\frac{1-e^z}{z-e^z}\,dz.$$ Moreover, $z=e^z$ implies: $$\|z\|=e^{\Re(z)}$$ hence all the complex solutions of $z=e^z$ lie on the curve $\gamma_+:(W(1),+\infty)\to\mathbb{C}$ $$ \gamma_+(r)=\log r+i\sqrt{r^2-\log^2 r}$$ or on the "conjugated" curve $$ \gamma_-(r)=\log r-i\sqrt{r^2-\log^2 r}.$$

Suppose we want to solve $z=e^z$. First let's start with Euler's Formula: $$ e^{x+iy}=e^x\cos(y)+ie^x\sin(y)\tag{1} $$ Thus, we want to solve $$ x=e^x\cos(y)\quad\text{and}\quad y=e^x\sin(y)\tag{2} $$

Locating the Roots

From $(2)$, we get $x^2+y^2=e^{2x}$. Therefore, the roots lie on the conjugate curves $$ y^2=e^{2x}-x^2\tag{3} $$ The roots of $z=e^z$ will lie on $(3)$ at the intersections with the curves $$ \cos(y)=xe^{-x}\tag{4} $$ Note that $(4)$ says $x\le0\implies\cos(y)\le0\implies y^2\ge\pi^2/4$. However, $e^{2x}-x^2\le1$ for $x\le0$. Therefore, all root have $x\gt0$.

Numerical Solution

To solve numerically, we can solve the equation $$ xe^{-x}=\cos\left(\sqrt{e^{2x}-x^2}\right)\tag{5} $$ for $x$ then use $(3)$ to get $y$.

As $x$ gets large, $xe^{-x}$ gets very small. This has two consequences: from $(3)$, we get $$ y=e^x\sqrt{1-\left(xe^{-x}\right)^2}=e^x+O\left(x^2e^{-x}\right)\tag{6} $$ and from $(4)$, we get that $y$ is very close to an odd multiple of $\pi/2$. Thus, for large $k$, we get $$ z_k\approx\log(k\pi+\pi/2)\pm i(k\pi+\pi/2)\tag{7} $$

Examples

Even for $k=0$, $(7)$ is not wildly wrong. $(7)$ gives $$ 0.45158271\pm1.57079633\,i $$ whereas using $(3)$ and $(5)$ to solve numerically gives $$ 0.31813151\pm1.33723570\,i $$ For $k=10$, $(7)$ approximates $$ 3.49610514\pm32.98672286\,i $$ while using $(3)$ and $(5)$ numerically gives $$ 3.49851521\pm32.88072148\,i $$ For $k=1000$, $(7)$ approximates $$ 8.05298504\pm3143.16344992\,i $$ while using $(3)$ and $(5)$ numerically gives $$ 8.05298751\pm3143.16088786\,i $$

For an answer within special functions, recall that the Lambert W function is defined as $W(z) e^{W(z)}=z$ (i.e. as the inverse function of $x e^x$.) Since the equation above can be written as $(-x)e^{-x}=-1$, we must have $x=-W(-1)$.

For consistency with Jack's result, note that the Lambert W function is an infinitely multi-valued function whose $n$-th branch is denoted by $W_n(x)$ for all integers $n$. WolframAlpha can be used to get a list of solutions for small $n$. Note that all the solutions are all complex and occur in complex conjugate pairs with $-W_{1-n}(-1)=-\overline{W_{n}(-1)}$.