I am supposed to find probability density of the random variable $Y$, which is an output of the following algorithm:
Gen $X \sim U(0,1)$
Gen $U \sim U(0,1)$
If $X<U$ then $Y:=X$ else $Y:=1-X$
return Y
In R:
X <- runif(10000, 0, 1)
U <- runif(10000, 0, 1)
res <- ifelse(X <= U, X, 1-X)
My attempts: $$ P(Y \le y) = \int \limits_0^1 P(Y\le y|U=u) f_U(u) du $$
$Y \in [0,1]$, so lets consider $y \in [0,1]$
$$ P(Y\le y|U=u) = P(Y\le y, \{ X \le u, X>u \}|U=u)= $$ $$ P(Y\le y, X \le u|U=u) + P(Y\le y, X>u |U=u) = $$ $$ P(X\le y, X \le u|U=u) + P(1-X\le y, X>u |U=u)= $$ $$ P(X\le y, X \le u|U=u) + P(X \ge 1-y, X>u |U=u)= $$ $$ P(X < \min(y,u)|U=u) + 1 - P(X \le 1-y, X<u |U=u)= $$ $$ P(X < \min(y,u)|U=u) + 1 - P(X < \min(1-y,u) |U=u)= $$ $$ \min(y,u) + 1 - \min(1-y,u) $$
$$ P(Y \le y) = \int \limits_0^1 \min(y,u) + 1 - \min(1-y,u) du = y - \frac{y^2}{2} +1 -\frac{1}{2} + \frac{y^2}{2} $$
$$ P(Y \le y) = y - 1/2 $$ That would mean that $Y \sim U[0,1]$ but if I perform the simulation in R it turns out that the PDF is equal to $-2x +2$
X <- runif(10000, 0, 1)
U <- runif(10000, 0, 1)
res<- ifelse(X <= U, X, 1-X)
hist(res, prob=TRUE, asp=1)
curve(-2*x+2, lwd=2, col="red", add=TRUE)
Where am I making a mistake?

One mistake -- Your transition from $$P(X\le y, X \le u|U=u) + P(X \ge 1-y, X>u |U=u)=$$ to $$P(X < \min(y,u)|U=u) + 1 - P(X \le 1-y, X<u |U=u)$$ is incorrect, because the complement of the event $(Z\ge a, Z>b)$ is not $(Z\le a, Z<b)$; it is $(Z<a \ \color{red}{\text{or}}\ Z\le b)$. A better way to proceed is to write $(Z\ge a, Z>b)=(Z>\max(a,b))$. From here you can continue with the complement device.
With this correction you will get the final form $$\min(y,u) + 1 - \max(1-y,u).$$ Integrating this over $u$ and then differentiating wrt $y$ will get you the correct form for the density.