Ok so I am very confused about the answer to this integral that I am getting. For some positive finite domain $A$ I want to solve
$$ \iint_{A} x^y dxdy $$
When I throw this in desmos' 3d plotter I have no issues and get finite values, but when I try to evaluate the integral by hand or with wolfram alpha I get the exact same solutions for the integral but when I try and evaluate the bounds this singularity at $x=1$ just "appears" and I can't explain why.
$$ \iint x^y dxdy =\int \frac{x^{y+1}}{y+1} dy=\int\frac{e^{(y+1)\ln(x)}}{y+1}dy $$ $$ =\text{Ei}((y+1)\ln(x)) $$
OR
$$ \iint x^y dxdy=\int \frac{x^y}{ln(y)} dx $$ $$ =\text{Li}(x^{(y+1)})=\text{Li}(e^{(y+1)\ln(x)})=\text{Ei}((y+1)\ln(x)) $$
Ok cool, but here is the problem, Ei(0) is a negative singularity which occurs when $x=1$ within $(y+1)\ln(x)$, so then if I tried to evaluate the integral of some region near it on my computer it diverges to negative infinity.
What's annoying is that my book, Table of Integrals, Series, and Products, agrees with me and so does wolfram alpha that the integral provided is in fact the solution. But say I want the square area from (0,0) to (1,1) it diverges or at least on my computer it ruins the computation in python. Furthermore, when I run this on desmos:
$$ z(x,y)=\int_0^x\int_0^y n^m dndm $$
It plots a function that clearly behaves as expected around 1. What am I messing up? Is there another analytical solution or am I confusing how to solve this right.
