var('x y')
f = x^y-y^x
max_x, max_y = 8, 8
range_x = (x, 0.01, max_x)
range_y = (y, 0.01, max_y)

p = plot([], dpi=150, gridlines=False)
p += region_plot(f>0, range_x, range_y, incol='lightblue', bordercol='gray', axes_labels=['$x$','$y$'])
e = exp(1)
p += text('$x^y > y^x$', [3.5, 5], rgbcolor='dimgray', fontsize='large')
p += text('$x^y > y^x$', [3.5, 1], rgbcolor='dimgray', fontsize='large')
p += text('$x^y < y^x$', [1.0, 4], rgbcolor='dimgray', fontsize='large')
p += text('$x^y < y^x$', [6.0, 4], rgbcolor='dimgray', fontsize='large')

point1 = (exp(1), pi)
point2 = (2, pi)
point3 = (exp(1), exp(1))

p += scatter_plot([point1], markersize=25, edgecolor='gray', facecolor='lightgray')
p += text(r'$(e,\,\pi)$', (point1[0], point1[1] + .25), fontsize='medium', rgbcolor='gray')

p += scatter_plot([point2], markersize=25, edgecolor='gray', facecolor='lightgray')
p += text(r'$(2,\,\pi)$', (point2[0] - .1, point2[1] + .25), fontsize='medium', rgbcolor='gray')

p += scatter_plot([point3], facecolor='red')
p += text(r'$(e,\,e)$', (point3[0], point3[1] - .5), fontsize='large', rgbcolor='black')

g_latex = r'$y = \exp\left(-W_0\left(\frac{-\ln x}{x}\right)\right)$'
g = lambda x: exp(-lambert_w(0, -ln(x)/x))
h_latex = r'$y = \exp\left(-W_{-1}\left(\frac{-\ln x}{x}\right)\right)$'
h = lambda x: exp(-lambert_w(-1, -ln(x)/x))
p += plot(g, range_x, thickness=2)
p += plot(h, (N(g(8)), 8), thickness=2, color='green')
p += text(g_latex, (6.2, 2.2), fontsize='x-large')
p += text(h_latex, (3.6, 7.0), fontsize='x-large', color='green')
p += text('$y = x$', (7.2, 6.4), fontsize='xx-large', color='black')
p += plot(x, range_x, linestyle=':', thickness=1, color='black')

show(p)
p.save('Plot_of_x^y_=_y^x.svg')