So what you want to do is evaluate the integral:
\int_0^\infty \int_0^\infty \lambda^2 e^{-\lambda(u+v)} du dv
by looking at the integrand in (x,y) space instead of (u,v) space, where x=u+v and y=uv. Is that correct?
If so, you're correct in the way you've approached the problem. You want to find the jacobian to transform du dv to dx dy, and then find the region in (x,y) space corresponding to the 1st quadrant in (u,v) space.
I think the problem is that x and y, as defined, do not uniquely determine a point (u,v). Because the defining equation are symmetric with respect to u and v, if (u,v) corresponds to the point (x,y), so does (v,u). In other words, x and y are bad coordinates.
That being said, this isn't a huge problem. As long as you restrict to one side of the line u=v (eg, the set of points (u,v) with 0 \leq u \leq v), x and y are still good coordinates. Luckily, the integral will be the same on each side of the line, so you can just compute the integral over one side and then double it.
The key point is to be careful with square roots. For example, x^2-4y = (u-v)^2 is true, but \sqrt{x^2-4y}= u-v is not true unless u>v, and if v>u this becomes \sqrt{x^2-4y}= v-u. Once you pick a side, you can pick the sign of the square root. In fact, by implicitly taking the positive root, you already integrated over the region u>v>0, and so your result of 1/2 makes perfect sense: you only integrated over half the first quadrant.