I read https://fse.studenttheses.ub.rug.nl/20174/1/bMATH_2019_DePooterJS.pdf
There is such an explanation in this article. By considering a $SL(2, \mathbb{Z})$action every geodesic on upper half plane can move the geodesic which is beginning point $u_{-\infty}$ and end points $u_{\infty}$ as follows: $0<\left|u_{-\infty}\right| \leq 1$ and $\left|u_{\infty}\right| \geq \infty$ and both are on opposite sides of the imaginary axis. Denote this set of geodesics by $A$. $A$ is the set of geodesic with points $u_{-\infty}$, $u_{\infty}$ belongs to $u_{-\infty} \in (0,1] \text { and } u_{\infty} \in (- \infty, -1]$ or $u_{-\infty} \in(-1,0] \text { and } u_{\infty} \in[1, \infty)$
I think that the power of $\tau=\left[\begin{array}{ll}1 & 1 \\ 0 & 1\end{array}\right]$ is not sufficient because geodesics with the distance between the starting point$u_{-\infty}$ and the endpoint $u_{\infty}$ is less than 2 cannot move to $A$.
How can we think about moving all geodesics in the upper half-plane to $A$ by considering approximate $SL(2, \mathbb{Z})$action?
I agree that using only $\tau$, it is not possible to make any two points on $\mathbb{R}$ lie in the respective intervals. This problem is possible exactly if it is possible to achieve that the points have distance at least one between each other.
If you also allow $$ \sigma = \begin{pmatrix}0 & -1 \\ 1 & 0\end{pmatrix}, $$ then it might be possible:
If the Euclidean distance between $u_{-\infty}$ and $u_{\infty}$ is at least 1, you are done. Otherwise apply $\tau$ until both points are inside the unit circle and then apply $\sigma$. This will increase their (Euclidean distance) as $|\sigma u_{-\infty} - \sigma u_{\infty} | = |u_{-\infty} - u_{\infty}|/{|u_{\infty}u_{-\infty}|} > |u_{-\infty} - u_{\infty}|$.
I work out the details for why this increases enough below:
We note that for $u \in (0,1)$ we have $\sigma(u) = -1/u$. We will also use the function $f= \tau \circ \sigma \circ \tau \circ \sigma$, which is described by $f(u) = 1/(1-b)$, (check this). In terms of distances, if $u,v \in (0,1)$, then $$ d(\sigma(u),\sigma(v)) = \frac{d(u,v)}{uv} $$ and $$ d(f(u),f(v)) = \frac{d(u,v)}{(1-u)(1-v)}. $$ We claim that no matter how $u,v$ are distributed on $(0,1)$, we can apply either $\sigma$ or $f$ to double their distance. Therefore by repeatedly applying $f$ and $\sigma$, we can increase the distance to above 1 at which point we can apply $\tau$ until the endpoints lie in the positions asked in the question. We show the claim by a case distinction:
Case 1: If $d(u,v) \geq 1/2$, then one of $u,v$ has to be at most $1/2$, so $uv \leq 1/2$. Then $d(\sigma(u),sigma(v)) \geq d(u,v)/(1/2) = 2d(u,v)$.
Case 2: $d(u,v) < 1/2$. Without loss of generality let $u \geq v$.
Case 2a: $u \geq 3/4$. Then $v > 1/4$ and $(1-u)(1-v) < 1/4 \cdot 3/4 = 3/8 < 1/2$. Then $d(f(u),f(v)) = d(u,v)/((1-u)(1-v)) > d(u,v)/(1/2) = 2 d(u,v) $.
Case 2b: $u < 3/4$. Then $v<1/4$ and $uv < 3/4 \cdot 1/4 < 1/2$. Then $d(\sigma(u),\sigma(v)) = d(u,v)/(uv) > 2 d(u,v)$.