Geodesics of hyperbolic plane I just found (or: I think that I found) the geodesics of the upper, closed half plane of $\mathbb R^2$. To verify my solution: Is it correct that the geodesics are the circles and lines which meet the unit circle at right angles? 
Thanks a lot for the confirmation (or the correction, if I'm wrong...)
Edit: Here is what I have so far (this should be correct, isn't it?):

 A: It is not correct that the geodesics of the upper half plane $\mathbb{H}=\{(x,y)\in\mathbb{R}^2:\;y>0\}$ with metric $g$ given by $g_{xx}=g_{yy}=1/y^2$ and $g_{xy}=g_{yx}=0$ the circles and lines which meet the unit circle at right angles. Instead, they are the circles and lines that meet the $x$-axis (the line $y=0$ in $\mathbb{R}^2$) at right angles.
A curve  $\gamma:I\subset\mathbb{R}\rightarrow \mathbb{H}:s\mapsto(x_1(s),x_2(s))$ is a geodesic if it satisfies the geodesic equation:
\begin{equation}
  \frac{d^2x_k}{ds^2} + 
  \sum_{i,j=1}^2\Gamma^k_{ij}\frac{dx_i}{ds}\frac{dx_j}{ds}=0,
\end{equation}
for $k=1,2$ and where $\Gamma^k_{ij}$ are the Christoffel symbols, which in this case are
\begin{equation}
\Gamma^1_{12}=\Gamma^1_{21}=\Gamma^2_{22}=-1/y;\;\;
\Gamma^2_{11}=1/y;\;\;\;
\Gamma^1_{11}=\Gamma^2_{12}=\Gamma^2_{21}=\Gamma^1_{22}=0.
\end{equation}
The geodesic equations are then the following system of differential equations (writing $x=x_1$ and $y=x_2$):
\begin{equation}
\frac{d^2x}{ds^2}-2\frac{1}{y}\frac{dx}{ds}\frac{dy}{ds}=0;\;\;\;
\frac{d^2y}{ds^2}+\frac{1}{y}
\left[\left(\frac{dx}{ds}\right)^2-\left(\frac{dy}{ds}\right)^2\right]=0
\end{equation}
For $\frac{dx}{ds}\neq 0$ the solutions satisfy: $x^2+y^2-ax=b$, where $a$ and $b$ are constants. Those are circles with center in the line $y=0$, and thus meeting that line at right angles. The solutions for $\frac{dx}{ds}=0$ are just vertical lines. 
A: Here I provide an alternative solution based on a physicist's convention.
In mathematical physics concerning differential geometry, e.g. general relativity, one often uses Killing vectors to help solving a geodesic. Each Killing vector $K$ gives a first integral (Einstein summation convention assumed)
$$
\frac{\mathrm{d} x^i}{\mathrm{d} \lambda} K_i = \text{const.},
$$
so that one can deal with a set of $n_K$ first-order differential equations, instead of the second-order and often non-linear equations, if the number of Killing vectors is $n_K$.
In your case, all the Killing vectors in Poincaré half-plane are given by
$$
T = \partial_1,\qquad
D = 2\left(x\,\partial_x+y\,\partial_y\right),\qquad
S = (y^2-x^2)\,\partial_x - 2xy\,\partial_y.
$$
Their musical isomorphisms $g(K,\cdot)$ are
$$
T^{\flat} = \frac{1}{y^2}\mathrm{d}x,\qquad
D^{\flat} = \frac{2}{y^2}\left(x\,\mathrm{d}x+y\,\mathrm{d}y\right),\qquad
S^{\flat} = \left(1-\frac{x^2}{y^2}\right)\,\mathrm{d}x - 2\frac{x}{y}\,\mathrm{d}y.
$$
One has therefore three first integrals
$$
t = \frac{\mathrm{d}x^i}{\mathrm{d}\lambda} T_i = \frac{1}{y^2} \frac{\mathrm{d}x}{\mathrm{d}\lambda},\\
d = \frac{\mathrm{d}x^i}{\mathrm{d}\lambda} D_i = \frac{2}{y^2} \left( x \frac{\mathrm{d}x}{\mathrm{d}\lambda} + y \frac{\mathrm{d}y}{\mathrm{d}\lambda}\right), \\
s = \frac{\mathrm{d}x^i}{\mathrm{d}\lambda} S_i = \left(1-\frac{x^2}{y^2}\right)\frac{\mathrm{d}x}{\mathrm{d}\lambda} - 2\frac{x}{y} \frac{\mathrm{d}y}{\mathrm{d}\lambda}.
$$
In this case, which is named as a maximally symmetric space in physics, the number of equations $n_K = N(N+1)/2 = 3$ is already larger than that of variables $N = 2$. Algebraic elimination gives
$$
(x^2 + y^2)t - x d - s = 0,
$$
which is essentially the same as coconut's answer.
An additional advantage of this method is, Killing vectors are realisation of symmetries, so that each integration constant now has a geometric meaning. Unfortunately as a physicist I am not able to express the meaning in mathematical terms.
A: It is possible to write down the solution in full parametric form. It is parametrized with the intial conditions: $x(0)=x_0$, $y(0)=y_0$, $x'(0)=u_0$,
$y'(0)=v_0$ and additionally the constraint $u_0^2+v_0^2=y_0^2$ is implied.
For the circle:
$x(\lambda) = x_0+\frac{u_0 y_0\sinh\lambda}{y_0\cosh\lambda-v_0\sinh\lambda}$
$y(\lambda)=\frac{y_0^2}{|y_0\cosh\lambda-v_0\sinh\lambda|}$
For $u_0=0$ it follows that $v_0=y_0$ and the straight line is recovered from the general solution:
$x(\lambda)=x_0$
$y(\lambda)=y_0 \exp(\lambda)$
