There is a geometric way of seeing that inversion works properly, which is based on considering not $1/z$ but $\bar{1/z}$. Then your $(x,y)$ gets mapped to $(\frac x{x^2+y^2},\frac y{x^2+y^2})$, which can be interpreted geometrically as follows. Let $O$ be the origin, and let $P$ be a point. Then $P$ gets sent to the unique point $P'$ on the ray $\vec{OP}$ such that $OP\cdot OP'=1$ (where $OP$ is the length of the segment joining $O$ and $P$). If $R$ is any point on the unit circle between $P$ and $P'$, then we can rewrite the equation as $\frac{OP}{OR}=\frac{OR}{OP'}$. Hence, triangles $OPR$ and $ORP'$ will be similar. You can then prove that triangles $OPQ$ and $OQ'P'$ are similar for any points $P$ and $Q$. (multiply to get $\frac{OP}{OR}\frac{OR}{OQ}=\frac{OQ'}{OR}\frac{OR}{OP'}$).
Once you have that it is easy to show that if a line $\ell$ does not go through the origin, it gets mapped to a circle that goes through the origin as follows. Suppose that $P$ is on $\ell$ such that $\vec{OP}$ is perpendicular to $\ell$. Then for any other point $Q$ on $\ell$ we have that triangle $OPQ$ has right angle $OPQ$, and hence triangle $OQ'P'$ has right angle $OQ'P'$ for all points $Q'$. But that implies that the points $Q'$ trace out a circle with diameter $OP'$. Hence, $\vec{OP}$, which was perpendicular to $\ell$ is perpendicular to the circle that $\ell$ got sent to. This shows full conformality since if you have two lines $\ell_1$ and $\ell_2$, you have three cases.
- $\ell_1$ and $\ell_2$ go through the origin -- they are fixed by inversion.
- $\ell_1$ goes thorugh the origin, but $\ell_2$ does not. Then $\ell_2'$ which is the perpenciular line to $\ell_2$ though the origin is fixed along with $\ell_1$ and perpendicular to the image-circle of $\ell_2$, hence the angle between $\ell_1$ and the circle equals the angle between $\ell_1$ and $\ell_2$.
- Both don't go through the origin, in which case take both of their perpendiculars $\ell_1$ and $\ell_2$.
There is also the analytic way of doing this.
If we consider $\mathbb C$ as the real plane $\mathbb R^2$, then what you're trying to show is that certain functions $f\colon \mathbb R^2\to\mathbb R^2$ are angle-preserving.
To be angle preserving at a particular point means that the derivative $df(z_0)$, which is a $2\times 2$ real matrix, preserves the angles between vectors. To see this, take two (smooth) curves $c_1\colon\mathbb R\to\mathbb R^2$ and $c_2\colon\mathbb R\to\mathbb R^2$. The angle between them at point $z_0=c(t_1)=c(t_2)$ is, as you said, the angle between the tangent lines, which is of course the angle between the tangent vectors $dc_1(t_1)$ and $dc_2(t_2)$ -- let's call them $v_1$ and $v_2$.
The function $f$ will map the two curves to the curves $f\circ c_1$ and $f\circ c_2$, whose tangent vectors will be $df(z_0)v_1$ and $df(z_0)v_2$. Thus, we want $df(z_0)$ to be a matrix $M$ that preserves the angles between vectors.
Fortunately, you're not asking how to classify all conformal mappings, so we just need to check that the derivatives of $z\to az+b$, $z\to 1/z$ and $z\to\bar z$ preserve angles at every $z_0$.
It turns out that the Cauchy-Riemann equations tell you that complex-differentiable functions (such as $z\to az+b$ and $z\to 1/z$) have derivates that look like $M=\left[\begin{matrix}a&b\\-b&a\\\end{matrix}\right]$ which since $\det M=a^2+b^2$ look exactly like $a^2+b^2\left[\begin{matrix}\cos(\theta)&\sin(\theta)\\-\sin(\theta)&\cos(\theta)\end{matrix}\right]$ for some angle $\theta$, which is exactly dilation by a factor of $a^2+b^2$ and rotation by angle $\theta$, which evidently preserve angles between vectors (for a quick proof, the rotation matrix is orthonormal and so actually preserves dot products, while dilation merely scales dot products by some constant).
To get the action of $\bar z$, note that it sends $x+iy$ to $x-iy$ and so its derivative is the matrix $\left[\begin{matrix}1&0\\0&-1\end{matrix}\right]$. This matrix does not preserve the dot product -- a short computation shows that it reverses the sign of the dot product. But that is ok since reversing the sign of the dot product is a reversal of orientation, i.e. the angle changes from $\theta$ to $180^o-\theta$, which for the purposes of being conformal is the same thing.
I kept making mistakes, so here is a pretty complete thing. Instead of guessing on the coefficients the way I intended, it is better to work with the cross ratio, pages 78-80. We start with the ordered triple $-1, 0, \frac{1}{2}$ and send these to $1,0,\infty$ by a Möbius transformation, namely $\frac{3z}{2z-1}.$ This takes the quadruple $-1, 0, \frac{1}{2},1$ to $1,0,\infty,3,$ so the cross ratio of the quadruple is $3.$
Next, we want a quadruple with some real $\alpha > 0$ given by $1,0, - \alpha, -1 - \alpha$ that also has cross ratio $3.$ That is, the Möbius trans given by $\frac{(1 + \alpha)z}{z + \alpha}$ must take $ -1 - \alpha$ to $3.$ So we get $\alpha = \sqrt 3 - 1$ and $1 + \alpha = \sqrt 3.$ So, $ \frac{\sqrt 3 \; z}{z + \sqrt 3 - 1} $ takes $1,0,1 - \sqrt 3, - \sqrt 3$ to $1,0,\infty,3.$ Which is good.
However, what we really want is to take $1,0,\infty,3$ to $1,0,1 - \sqrt 3, - \sqrt 3,$ which requires the inverse item, namely $\frac{(\sqrt 3 - 1) \; z}{-z + \sqrt 3}.$
Taking the composition and dividing all entries by $\sqrt 3,$ we get
$$ f(z) = \frac{(3 - \sqrt 3) \; z \; }{(2 - \sqrt 3) \; z - 1} $$
This gives us $$f(-1) = 1, \; \; f(0) = 0, \; \; f(\frac{1}{2}) = 1 - \sqrt 3, \; \; f(1) = - \sqrt 3.$$
Let's see, the determinant is negative, real entries, so it does fix the real axis but takes positive imaginary part to negative imaginary part.
The common center is the real number $$ \frac{1 - \sqrt 3}{2}. $$ If you wish, you may subtract that off to move the center(s) to the origin, then multiply by a real number to adjust the radii. Thus it is possible to get the unit disk back to where it was.
Best Answer
You may assume that all four circles are centered at the origin. The $x$-axis $l_1$ and the $y$-axis $l_2$ can be considered as circles that intersect the first pair of circles at $90^\circ$. They will be mapped by $f$ onto two circles that intersect the second pair of circles at $90^\circ$, and it is easy to see that this is only possible if the circles $f(l_i)$ are again lines through the origin. As $l_1$ and $l_2$, as well as $f(l_1)$ and $f(l_2)$, intersect at $0$ and $\infty$ it follows that $f$ either keeps $0$ and $\infty$ fixed or interchanges these two points. This in turn implies $f(z)=c z$ or $f(z)=c/z$ for a suitable $c\ne 0$. In both cases the ratio between the larger and the smaller radius of the two circles stays the same.