As far as I understand it the closing lemma implies that closed geodesics on surfaces of negative curvature are dense. So: how can they be constructed in general?
A concrete answer that dovetails with the construction of such surfaces with constant negative curvature and genus $g$ from regular hyperbolic $(8g-4)$-gons along lines indicated by Adler and Flatto and gives the endpoints of the geodesics in the Poincaré disk model would be ideal. More useful still would be a way to construct all the closed geodesics that cross the boundaries of translates of the fundamental $(8g-4)$-gon some specified number of times (I am pretty sure this ought to be a finite set, but I couldn't say why off the top of my head).
Best Answer
If you think of your surface as the upper half plane modulo a group of Moebius transformations $G$, start by representing each of your Moebius transformations $ z \longmapsto \frac{az+b}{cz+d}$ by a Matrix.
$$A = \pmatrix{ a & b \\\ c & d}$$
And since only the representative in $PGL_2(\mathbb R)$ matters, people usually normalize to have $Det(A) = \pm 1$.
The standard classification of Moebius transformations as elliptic / parabolic / hyperbolic (loxodromic) is in terms of the determinant and trace squared. You're hyperbolic if and only if the trace squared is larger than $4$. Hyperbolic transformations are the ones with no fixed points in the interior of the Poincare disc, and two fixed points on the boundary, and they are rather explicitly "translation along a geodesic".
Elliptic transformations fix a point in the interior of the disc so they can't be covering transformations. Parabolics you only get as covering transformations if the surface is non-compact, because parabolics have one fixed point and its on the boundary -- if you had such a covering transformation it would tell you your surface has non-trivial closed curves such that the length functional has no lower bound in its homotopy class.
So your covering tranformations are only hyperbolic. That happens only when $tr(A)^2 > 4$. So how do you find your axis? It's the geodesic between the two fixed points on the boundary, so you're looking for solutions to the equation:
$$ t = \frac{at+b}{ct+d}$$
for $t$ real, this is a quadratic equation in the real variable $t$. If I remember the quadratic equation those two points are:
$$ \frac{tr(A) \pm \sqrt{tr(A)^2 - 4Det(A)}}{2c}$$
Is this what you're after?