Geodesics of hyperbolic plane

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.

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.