The error stems from the fact that you use $ux uy$ instead of just $uy$, which introduces a non-random $x$-dependent scaling factor for the uniform random $uy$ when sampling for $y$ and also the use of $\text{cdf}_{xy}$ instead of $\text{cdf}_{y|x}$. The remaining portion seems to be correct.
Inverse Transform Sampling
If $u \sim U(0,1)$, and $x \sim P_x$, then
$$\text{cdf}_x=\int_0^xP_xdx$$
$$\text{if}\quad \bar x = \text{cdf}_x^{-1}u\quad \text{then}\quad \bar x \sim P_x$$
The Ideal way
(Let $u_1$ and $u_2$ be the uniform random variables and $\text{cdf}$ be the cumulative distribution function.)
$$\text{cdf}_{xy}(x,1) = u_1$$
$$x =\text{cdf}^{-1}u_1 =6-\sqrt{36-11u_1}\quad\implies \quad x\sim \left(P_{x}=\int_0^1P_{xy}dy\right)$$
$$\text{cdf}_{y|x}(x,y) = u_2\quad \text{for} \quad x=6-\sqrt{36-11u_1}\quad\implies\quad \text{for given }x,\quad y\sim \left(P_{y|x}=\frac{P_{xy}}{P_x}\right)$$
$$\text{cdf}_{y|x} = \int_0^y P_{y|x} dy = \int_0^y \frac{P_{xy}}{P_x} dy$$
$$\therefore x,y \sim (P_{y|x}P_x=P_{xy})$$
The Incorrect Analysis
The initial part in your analysis is correct, and hence $x$ follows the marginal distribution it is supposed to follow, which is apparent from the fact that the values for $\text{cdf}(x,1)$ matches the expected values. The second part is faulty due to the $u_1$ term as follows:
$$\text{cdf}_{xy}(x,y) = u_1u_2\quad\text{for}\quad u_1=\frac{1}{11}(12x-x^2)=f(x)$$
$$\therefore \text{for given }x,\quad \frac{\text{cdf}_{xy}(x,y)}{f(x)} = u_2$$
So, because of the $f(x)$, we end up using a different $\overline{\text{cdf}}_{y|x}$ to sample $y$ instead of the intended cumulative density function given a particular $x$. To get the $\overline{\text{cdf}}_{y|x}$, we assume this to be coming from a faulty $\overline{P}_{y|x}$:
$$\overline{P}_{y|x} = \frac{\partial}{\partial y}(\overline{\text{cdf}}_{y|x})$$
$$\overline{P}_{y|x} = \frac{\partial}{\partial y}\frac{8xy+2x^2y+4xy^2-3x^2y^2}{12x-x^2} = \frac{-6xy+2x+8y+8}{12-x}$$
$$\therefore \overline P_{xy} = \overline P_{y|x}P_x$$
$$\overline{\text{cdf}}_{xy} = \int_0^y \int_0^x \overline P_{xy} dx dy$$
This, as solved in WolframAlpha is a very complicated function, but for our purposes, we can just take the $x$ integral from 0 to 1.
$$\overline{\text{cdf}}_{xy}(1,y) = -\frac{2}{11}y\left(\frac{31y}{2}-384(y-1)\coth ^{-1}(23)-21\right)$$
Evaluating this at different values of $y$ will give you the desired faulty values that you are getting in the sampling case above.
Addendum
I did a slight naming error in the Ideal way above, and I have now fixed it. Instead of taking $\text{cdf}_{xy}$ we needed $\text{cdf}_{y|x}$, which is what I assumed from that step onwards but wrote incorrectly. Corrected. That said,
Consider we have chosen u1 such that x=0.5. Then, cdf(x,1) is 0.522.
So, for any value 0≤y≤1, cdf(x,y) is smaller or equal to 0.522. But
what if we sampled u2 greater than 0.522? Then, y would need to be
greater than 1, which is invalid. That's why I thought that
normalizing with cdf(x,1) would be necessary. Or, alternatively
sampling u2 from [0,cdf(x,1))=[0,u1)
If, $u_1 \approx 0.522$, that would give us $x \approx 0.5$ from the marginal $x$ sampling (the $\text{cdf}_{xy}(x,1) = \text{cdf}_x = \int_0^x P_x dx$). That means, for $y$, we would set $\text{cdf}_{y|x}(0.5,y) = u_2$. The substituted $0.5$ already takes care of the normalization such that at $u_2=1$, we just need $y=1$
In our case, $\text{cdf}_{y|x} = (2y+xy+y^2-1.5xy^2)/(3-0.5x)$. Substitute $x=0.5$, to get $(2.5y+0.25y^2)/(2.75)=1$, which gives $y\in\{1,-11\}$.
Best Answer
The reasoning seems basically correct. E.g., in the scalar case, given some invertible $g:x\mapsto y$ and a lebesgue density $f_X(x)$ for $x$, a density for $y$ is some function $f_Y$ on the image of $g$ such that $\int_A f(x)dx=\int_{f(A)} f_Y(y)dy$, i.e., $f_Y(y)=dx/dy=1/g'(x)=1/g'(g^{-1}(y))$. For $x\in \mathbb{R}^N$ and $g$ invertible, we get $f_y(y)=|g'(g^{-1}(y))|.$ So if $g$ is linear, then $f_Y$ is constant.
If $g$ is a singular linear transformation $g(X)$ need not be uniform, e.g., let $X$ be uniform on a triangle in a plane and then the projection onto an axis is non-uniform. That is why the author makes the assumption about the points being in general position.