This is just a partial answer.
Since I faced this specific problem a few years ago for $W_0(x)$, what I found as "best" are Padé approximants.
The simplest ones are
For $\color{red}{-\frac 1e \leq x \leq -\frac 1{2e}}$
$$W_0(x) \approx \frac{-1+\frac{14\sqrt{2}}{45} \sqrt{e x+1}+\frac{301}{540} (e x+1)}{1+\frac{31\sqrt{2}}{45} \sqrt{e x+1}+\frac{83}{540}
(e x+1)}\tag 1$$
For $\color{red}{-\frac 1{2e} \leq x \leq 0}$
$$W_0(x) \approx \frac{x+\frac{4 }{3}x^2}{1+\frac{7 }{3}x+\frac{5 }{6}x^2}\tag 2$$
For sure, for more accuracy, I built similar expressions with more terms.
Edit
For the other branch, you can use
$$W_{-1}(x)=L_1-L_2+\frac{L_2}{L_1}+\frac{L_2(-2+L_2)}{2L_1^2}+\frac{L_2(6-9L_2+2L_2^2)}{6L_1^3}+\frac{L_2(-12+36L_2-22L_2^2+3L_2^3)}{12L_1^4}+\frac{L_2(60-300L_2+350L_2^2-125 L_2^3+12 L_2^4)}{60L_1^5}+\cdots$$ where $L_1=\log(-x)$ and $L_2=\log(-L1)$
Update
Calling $f(x)$ and $g(x)$ the expressions given in $(1)$ and $(2)$, let us consider the error function
$$\Phi(a)=\int_{-\frac 1 e}^a \left(W(x)-f(x)\right)^2+\int^0_a \left(W(x)-g(x)\right)^2$$ and get the following values
$$\left(
\begin{array}{cc}
a & 10 ^{10} \,\Phi(a) \\
-0.350 & 472433 \\
-0.325 & 72500 \\
-0.300 & 13455 \\
-0.275 & 2646 \\
-0.250 & 523 \\
-0.225 & 109 \\
\color{red}{ -0.200} &\color{red}{ 44} \\
-0.175 & 57 \\
-0.150 & 100 \\
-0.125 & 168 \\
-0.100 & 269 \\
-0.075 & 411 \\
-0.050 & 605 \\
-0.025 & 863 \\
0.000 & 1197
\end{array}
\right)$$
Consider that you look for the zero's of function
$$f(x)=x^x-x-1$$
Its first derivative $f'(x)=x^x (\log (x)+1)-1$ cancels at $x=1$ and the second derivative test $f''(1)=2$ shows that this is a minimum.
Build a Taylor expansion to get
$$f(x)=-1+(x-1)^2+\frac{1}{2} (x-1)^3+\frac{1}{3} (x-1)^4+O\left((x-1)^5\right)$$ Using series reversion, then
$$x=1+\sqrt{y+1}-\frac{y+1}{4}-\frac{1}{96} (y+1)^{3/2}+O\left((y+1)^2\right)$$ where $y=f(x)$. Making $y=0$, this gives as an approximation
$$x=\frac{167}{96}\approx 1.73958 $$ To polish the root, use Newton method starting with this estimate. The iterates will be
$$\left(
\begin{array}{cc}
n & x_n \\
0 & 1.739583333 \\
1 & 1.778584328 \\
2 & 1.776779132 \\
3 & 1.776775040
\end{array}
\right)$$
Edit
If we make the first expansion $O\left((x-1)^n\right)$ and repeat the inversion series, we generate the sequence
$$\left\{2,\frac{7}{4},\frac{167}{96},\frac{175}{96},\frac{160
379}{92160},\frac{3687}{2048},\frac{12144341}{6881280},\frac{110221693}{61931520
},\frac{211659504277}{118908518400}\right\}$$
We can also use $x_0=2$ and use high order iterative methods. For order $4$, that is to say one level after Householder method, we have
$$x=2\,\frac {4575+67460 a+299400 a^2+558920 a^3+463660 a^4+141128 a^5} {6655+86720 a+352260 a^2+615000 a^3+483960 a^4+141128 a^5 }$$ where $a=\log(2)$.
This gives, as another approximation, $x=1.776779506$.
Best Answer
Let's consider solutions in closed form.
$$e^x\ln(x)=a$$
1.) No algebraic solutions except $1$ for algebraic $a$
According to Schanuel's conjecture, we have $\forall\ 0,1\neq x\in\overline{\mathbb{Q}}\colon$ $\text{trdeg}_\mathbb{Q}\mathbb{Q}\left(x,\ln(x),e^x\right)=2$. That means: $\ln(a_0),e^{a_0}$ are $\overline{\mathbb{Q}}$-algebraically independent for each algebraic $a_0\neq 0,1$. That means, your equation cannot have algebraic solutions except $1$ for algebraic $a$.
2.) No general elementary solutions
According to the theorem of [Ritt 1925], that is proved also in [Risch 1979], the function $f\colon x\mapsto e^x\ln(x)$ seems not to have partial inverses that are elementary functions for non-discrete domains of $f$.
But that doesn't say if solutions in the elementary numbers are possible for single $a$.
[Risch 1979] Risch, R. H.: Algebraic Properties of the Elementary Functions of Analysis. Amer. J. Math. 101 (1979) (4) 743-759
[Ritt 1925] Ritt, J. F.: Elementary functions and their inverses. Trans. Amer. Math. Soc. 27 (1925) (1) 68-90
3.) Solutions in terms of Hyper Lambert W
$$e^x\ln(x)=a$$
$$x\to e^z: e^{e^z}(z+2k\pi i)=a\ \ \ \ \ (k\in\mathbb{Z})$$
$$-\pi<\text{Im}(z)\le\pi:$$
$$e^{e^z}z=a$$
This equation can be solved by Hyper Lambert W:
Galidakis, I. N.: On solving the p-th complex auxiliary equation $f^{(p)}(z)=z$. Complex Variables 50 (2005) (13) 977-997
Galidakis, I. N.: On some applications of the generalized hyper-Lambert functions. Complex Variables and Elliptic Equations 52 (2007) (12) 1101-1119