There are two problems involved here: One, is a general local root finding problem, for which many method's (such as Newton's) can be used. The second, is the problem of finding all (at most three) solutions.
One method I could suggest is the "Homotopy Continuation Method". To use it, first find all three solutions for some $a> e^e\sim15.15$. Now, vary $a$ slightly $a\to a+\Delta a$, and use the previous roots $x_i$ as initial values for your solver. Continue this process literately, until you find the solution curves for all $a$.
You can see how the solutions converge around the bifurcation $a= e^e$, as calculated using the attached python script.
The bifurcation is at $a= e^e$, since three solution can only exist when $df/dx=0$. This happens when $a^x = \log^2 (a) x$, and since the Lambert W function is not real valued below $-1/e$, we get that $a\ge e^e$.
from math import log
return a**(-x) + log(x) / log(a)
return x * (log(a) + (a**x) * log(x) ) / ( (a**x) - x * log(a) * log(a))
# find initial roots for a = 16
a = 16.0
roots = 
for x in [x * 0.1 for x in range(1, 10)]:
y = x
while (y > 0 and abs(f(y,a)) > 1.0e-14):
y = y - f_over_fprime(y,a)
if (y < 0): y = 0
found = False
for x in roots:
if (x == round(y,14)): found = True
if (found == False): roots.extend([round(y,12)])
print "Roots for a = 16: ", roots
for a in [a * 0.1 for a in range(160, 20, -1)]:
y = roots
for i in [0,1,2]:
while (y[i] > 0 and abs(f(y[i],a)) > 1.0e-14):
y[i] = y[i] - f_over_fprime(y[i],a)
if (y[i] < 0): y[i] = 0
roots = y
print "Roots for a = ", a, "are: ", roots