Part 1
Let's solve this from first principles and then reflect on what the solution might teach us.
First Step: Finding the Marginal and Joint Density Functions
The marginal density function of $X$ is
$$f_X(x) = I_{[0,1]}(x),$$
the indicator function for its interval of support $[0,1]$.
The conditional density of $Y$ is
$$f_{Y|x}(y) = \phi(y;x,|x|)= \phi(y;x,x)$$
where $\phi(x;\mu,\sigma)$ is the Normal PDF with mean $\mu$ and standard deviation $\sigma$ (and $|x|=x$ since $0\le x\le 1$). Therefore the joint density function is
$$f_{X,Y}(x,y) = I_{[0,1]}(x) \phi(y;x,x).$$
The marginal density of $Y$ is obtained by integrating out $x$,
$$f_Y(y) = \int f_{X,Y}(x,y) dx = \int_0^1 \phi(y;x,x) dx.$$
On the principle of Just-in-time-computation, I will not evaluate these integrals until that might be needed. (It turns out it is not needed at all!) Note, therefore, that no calculations have yet been performed: everything is just applying definitions.
Warm-up: Reproducing the Moments in the Question
Let's proceed to compute moments. The moments of $X$ are the moments of the uniform distribution, which are easily calculated or looked up. (The mean is $\mu_X^{(1)}=1/2$ and the mean square is $\mu_X^{(2)}=1/3$, for instance.) The $k^\text{th}$ moment of $Y$, $\mu_Y^{(k)}$, can be computed as
$$\mu_Y^{(k)} = \int y^k f_Y(y) dy = \int y^k \int_0^1 \phi(y;x,x) dx dy.$$
Let's interchange the order of integration (Fubini's Theorem) to obtain
$$\mu_Y^{(k)} = \int_0^1 \left(\int y^k \phi(y;x,x) dy\right) dx.$$
The inner integral computes the $k^\text{th}$ moment of a Normal variable $Y$ with mean $x$ and variance $x^2$, so it is immediate that this moment is $x$ when $k=1$ and $x^2 + (x)^2 = 2x^2$ when $k=2$. Whence
$$\mu_Y^{(1)} = \int_0^1 (x) dx = \frac{1}{2}$$
and
$$\mu_Y^{(2)} = \int_0^1 (2x^2) dx = \frac{2}{3}.$$
We conclude the variance of $Y$ is $\mu_Y^{(2)} - \left(\mu_Y^{(1)}\right)^2 = 2/3 - (1/2)^2 = 5/12$.
Solution
Having established that this method of direct integration, using Fubini's Theorem, reproduces the calculations in the question, let's apply it to the covariance:
$$\mathbb{E}(XY) = \int xy f_Y(y) dy = \int_0^1 x \left( \int y \phi(y;x,x) dy \right) dx = \int_0^1 x(x)dx = \frac{1}{3}.$$
Therefore
$$\text{Cov}(X,Y) = \mathbb{E}(XY) - \mathbb{E}(X)\mathbb{E}(Y) = \frac{1}{3} - \left(\frac{1}{2}\right) \left(\frac{1}{2}\right)= \frac{1}{12}.$$
Reflection
It should now be evident that the reversal of the order of integration is really the same thing as the Law of Iterated Expectation (LIE). The strategy was to compute the covariance in terms of the raw moments $\mathbb{E}(XY)$, $\mathbb{E}(X)$, and $\mathbb{E}(Y)$ by exploiting the knowledge of the conditional moments of $Y$, given $X=x$, and integrating those over $x$. Thus we could have invoked the LIE at the outset in the form
$$\mathbb{E}(XY) = \mathbb{E}(\mathbb{E}(XY | X)) = \mathbb{E}(X\,\mathbb{E}(Y | X)) = \mathbb{E}(X(X)) = \mathbb{E}(X^2) = \frac{1}{3}$$
with a similar (but easier) manipulation to find $\mathbb{E}(Y)$. Comparison with the previous integral expression shows that the crucial idea is to factor $X$ out of the conditional expectation:
$$\mathbb{E}(XY | X) = X\,\mathbb{E}(Y | X).$$
That may be the main thing worth remembering from this exercise.
Although these results can be used (in a straightforward manner) to derive an iterated law for covariances, comparable to the law for variances quoted in the question, I find such rules to be a burden on the memory. By identifying and focusing on fundamental principles and ideas, we only need to remember and learn a few things--and we can learn them well. Everything else can then easily be derived as needed.
Part 2
Let $xz = y$, entailing $dy = x dz + z dx$. Since
$$\eqalign{
\phi(y;x,x)dxdy &= \frac{1}{x}\phi\left(\frac{y}{x};1,1\right) dxdy \\
&= \frac{1}{x}\phi\left(z;1,1\right) dx(x dz + z dx) \\
&= x\frac{1}{x}\phi\left(z;1,1\right) dx dz \\
&= \phi\left(z;1,1\right) dx dz,
}$$
the joint distribution easily factors as
$$f_{X,Y}(x,y)dxdy = I_{[0,1]}(x) \left[\phi(z;1,1)dxdz\right] = \left[I_{[0,1]}(x)dx\right]\, \left[\phi(z;1,1)dz\right].$$
Since all probability calculations can now be carried out separately for $x$ and $z$, it follows that $X$ and $Z=Y/X$ are independent, QED.
(The manipulations with differentials are interpreted as wedge products. This is fully explained at https://stats.stackexchange.com/a/154298.)
As a check, we may simulate a large number of independent realizations of $(X,Y)$ and plot them. The following carries this out with R:
n <- 1e4
x <- runif(n)
y <- rnorm(n, x, x)
plot(x, y/x, cex=1/3, col="#00000020")
In any vertical slice the distributions of the $y_i/x_i$ look approximately the same. (They are all standard Normal with mean $1$, as we can see from their PDF $\phi(z-1)$.) That is what independence of $Z$ and $X$ means: the distribution of $Z$ does not vary with the value of $X$ (nor does the distribution $X$ vary with the value of $Z$).
