Olá, Bruno ...o autor fez o gráfico em polares, com o ParametricPlot3D
ou similar, daí o Mesh já vem com uma opção default que deixa tudo cheio
de radiais (theta constante) e circunferências levantadas para o
gráfico, note que temos
z = x y /(x^2 + y^2)^(1/2)
...em polares teremos
x = r costheta,
y = r sentheta,
z = simplificando = r costheta sentheta,
daí, para ficar igual ao livro, entramos com a parametrização
ParametricPlot3D[{r Cos[theta], r Sin[theta], r Cos[theta] Sin[theta]},
{r, 0, 1}, {theta, 0, 2pi} etc.]
e você não precisa do RegionFunction e outros controles de domínio de
plotagem, como estava usando, a limitação do r, nas polares, já implica
em região de plotagem cilíndrica.
Executei no seu note, que virou 20200423c.nb, logo após suas execuções o
comando:
ParametricPlot3D[{{r Cos[theta], r Sin[theta], r Cos[theta] Sin[theta]},{r
Cos[theta], r Sin[theta], -1}}, {r, 0, 1}, {theta, 0, 2pi} etc.],
pois é, este comando aceita duas parametrizações entre chaves (quase
tudo no mathematica aceita lista de coisas entre chaves) ...e perceba,
na segunda parametrização coloquei a coordenada z fixa e igual a -1,
para aparecer a sombra do gráfico, com o mesmo Mesh, no plano z = -1, lá
embaixo ...acho que fica bonitinho, é uma representação do que seria um
domínio circular da função com as divisões típicas das polares, o
gráfico acima deste domínio pode ser entendido como o domínio que foi
levantado para o espaço e entortado, pela alteração da coordenada z
...legal, né?
Saudações. Márcio.