I think you can understand why LASSO selects at most $n$-features intuitively if you think about the contours of $rss = |X\beta - y|_2^2$ and $l_1 = |\beta|_1$ on the parameter space $\mathbb R^p$. You can think about finding the minimal $rss + l_1$ solution as looking at the 2 contours of larger and larger values for $rss$ and $l1$$l_1$ until they intersect for the first time. Their growth with respect to each other needs a specific ratio, but that's irrelevant for the intuition, in fact in your head you might want to keep the $l_1$-contour constant.
The $l1$$l_1$-contour will be a generalized octahedron (https://en.wikipedia.org/wiki/Octahedron) with the corners on the axis. For your example with $p = 2$ that's just a square with sides 45 degrees rotated against the axis. For $rss$ let's first think about $rss = 0$, so the solution space to $X\beta = y$. That's a $p-n$-dimensional affine subspace. In your example that's the line $\beta_2 = -0.5\beta_1 + 2.5$. $rss$ is just distance squared so constant $rrs \implies$ constant distance. Therefor the contours will be parallel(or a elliptic, symmetric arrangement of parallel spaces) to the original subspace.
By going through one of the corners, which is on an axis, so the other parameter will be $0$.
By being parallel to one of the sides and touching the whole side, in which case you have non unique solutions (https://www.stat.cmu.edu/~ryantibs/papers/lassounique.pdf) This is partially discussed in some of the threads you linked. If you change $x2 = 1$$x_2 = 1$ in your example, then you will have an obvious case.
Generalizing: the $l1$$l_1$- contourcontour is made up of $d$-dimensional edges ($d = 0$ is a corner, $d = 1$ is an edge, $d = 2$ is a 2d-surface like the face for $p = 3$, ...) and hitting that $d$-dimensional edge correspond to selecting $d+1$ features. The $rss$-contour covers $p-n$-dimension and the $n$-dimensional edge, that would select $n+1$ features, covers the remaining dimensions, so there is no space to not be parallel.