I am stucking at the step which is determining new boundaries for variables after changing coordinate system. The question is : If $0\le x\le 1$ and $0\le y \le 1$ then after substituting $s=x+y$ and $p=xy$ then what is the new bound for $s,p$.
My attempt is $0\le x\le 1$ and $0\le y \le 1$, so $0\le p\le 1$. And since for existence of $x,y$ then $s^2\ge 4p$ so $2\sqrt{p}\le s$. But I don't know the upper bound for $s$.
Can I ask for help? Thank you so much.
Notice that when you compute the inverse Jacobian
$$J^{-1} = \begin{vmatrix}1 & 1 \\ y & x\end{vmatrix} = |x-y| \implies J = \frac{1}{|x-y|}$$
In other words, the Jacobian for the transformation has a big honking singularity that divides your square into two - because on both sides of the diagonal, the square maps to two copies of the same region.
How to deal with this? Well, first off the coordinate transformation you propose is invalid for your region - you would have to split it up into two integrals:
$$\int_0^1\int_0^1 f(x+y,xy)dxdy = \int_0^1\int_0^xf(x+y,xy)dydx + \int_0^1\int_0^yf(x+y,xy)dxdy$$
Let's focus only on the first integral for now which is defined by the triangle with the boundaries $y=0$, $x=1$, and $y=x$. Using the coordinate transformation, each of these boundaries becomes
$$\begin{cases}y = 0 \implies p = 0 \\ x = 1 \implies p = s-1 \\ y=x \implies 4p = s^2\end{cases}$$
which when drawn out looks like the following region in the $sp$ plane:
Verify for yourself that the other triangle also maps to this same region. To keep the number of integrals as minimal as possible, it's obvious the $s$ integral has to be done first (but this may not be possible depending on what your integral is, you can always do $p$ first and then break these two integrals down into four).
All that's left is to compute the Jacobian in $sp$ coordinates. We secretly already did most of this when we compute the boundary transformation for $y=x$
$$s^2-4p = (x-y)^2 \implies \sqrt{s^2-4p} = |x-y|$$
with no ambiguity in sign because on our region of integration in the $sp$ plane, $s^2$ is always larger than $4p$. This means our final transformation is
$$\int_0^1 \int_{2\sqrt{p}}^{p+1} f(s,p)\frac{dsdp}{\sqrt{s^2-4p}} + \int_0^1 \int_{2\sqrt{p}}^{p+1} f(s,p)\frac{dsdp}{\sqrt{s^2-4p}} = \int_0^1 \int_{2\sqrt{p}}^{p+1} f(s,p)\frac{2dsdp}{\sqrt{s^2-4p}}$$
In other words, we get a factor of $2$ in our Jacobian. We would have never discovered this factor of $2$, nor would we have found the third boundary of the region $s^2 = 4p$ had we not split the region up properly at $x=y$ in the first place.