Plotting geodesics in $SL_2(\mathbb{Z})\backslash\mathbb{H}$
In here there are functions that allow to plot geodesics in the upper half plane (namely, half circles perpendicular to the x axis and vertical lines). Are there similar functions that do the same thing but only in the fundamental domain of $SL_2(\mathbb{Z})$, so that the geodesic is folded up into the domain $|x|<\frac{1}{2}$ and $x^2+y^2>1$?
What I would like to have is something like the picture in page 5 of this paper.
There is no straightforward function to do it. However a simple algorithm is not hard to implement. You just need to look at the next intersection with the boundary, and apply the appropriate isometry... though the intersection feature is broken: see #23427.