I'm well aware that the Möbius transformations that take the unit disk to itself, $f:D \to D$, are given by $$ f(z) = \frac{e^{i \theta}(z-\alpha)}{1-\bar{\alpha}z}, $$ where $\theta\in [0,2\pi)$ and $|\alpha|<1$. Moreover, it's easy to see that regardless of the value for $\theta$, we have $f(\alpha) = 0$.
I'm also aware that Möbius transformations, in general, are completely determined by where they map three points; and that if you want to find the particular transformation that sends $z_1 \mapsto w_1$, $z_2\mapsto w_2$, and $z_3\mapsto w_3$, it is given by $$ f(z) = \frac{az+b}{cz+d}, $$ where $$ a = \det \begin{pmatrix} z_1 w_1 & w_1 & 1\\ z_2 w_2 & w_2 & 1\\ z_3 w_3 & w_3 & 1 \end{pmatrix}, \quad b = \det \begin{pmatrix} z_1 w_1 & z_1 & w_1\\ z_2 w_2 & z_2 & w_2\\ z_3 w_3 & z_3 & w_3 \end{pmatrix} $$ $$ c = \det \begin{pmatrix} z_1 & w_1 & 1\\ z_2 & w_2 & 1\\ z_3 & w_3 & 1 \end{pmatrix}, \quad d = \det \begin{pmatrix} z_1 w_1 & z_1 & 1\\ z_2 w_2 & z_2 & 1\\ z_3 w_3 & z_3 & 1 \end{pmatrix}. $$
So naturally, shouldn't there be some way of explicitly constructing a Möbius transformation that preserves the unit disk and takes $z_i \mapsto w_i$, for $i=1,2,3$, with $|z_i|,|w_i|<1$ given?
The problem that I'm running into is that the formula given above in terms of determinants is not preserving the unit disk. I'm using Mathematica to plot the image of a number of coordinate curves, and it's clear that the disk is not being preserved.
What am I missing? How should I choose my values of $\alpha$ and $\theta$ to make the correct transformation? In this case there are two degrees of freedom, yet in the other case there seem to be four. How do I reconcile this? Trying to equate the two formulas, and setting up a system of equations like $$ a = e^{i\theta}, \quad b = -\alpha e^{i\theta}, \quad c= -\bar{\alpha}, \quad d=1 $$ doesn't seem very helpful either.
What you're asking for is impossible.
First of all, we know that some Möbius transformations do not take the unit disc to itself, for example $f(z)=\frac{z+1}{0z+1} = z+1$.
On the other hand, if we pick any three distinct points $z_1,z_2,z_3$, and then compute $w_i=f(z_i)$ for $i=1,2,3$, it follows that $f$ is the unique Möbius transformation such that $w_i=f(z_i)$ for $i=1,2,3$.
So, using $f(z)=z+1$, let $z_1=1$, $z_2=2$, $z_3=3$. Then compute $w_1=2$, $w_2=3$, $w_3=4$. With those values of the $z_i$'s and $w_i$'s, what you want is doomed to failure.
ADDED LATER:
Regarding your additional question asked in the comments regarding the Appolonian gasket, what that link is saying is that in this image, once you have chosen the outer circle $C$ and three points $x,y,z$ on the circle $C$, the next three circles inside $C$, tangent to each other, and tangent to $C$ at $x,y,z$, are determined by the choices of $C,x,y,z$. Also, all subsequent circles in the construction of the gasket are also determined. Then comes this sentence:
What this means is that if instead you had instead chosen three other points $x',y',z'$ on the circle $C$, then the Möbius transformation that takes $x,y,z$ to $x',y',z'$ must also take $C$ to $C$. And, it must take the first three circles of the $x,y,z$ gasket to the first three circles of the $x',y',z'$ gasket. And this continues, by induction, as you go down to deeper and deeper levels of circles.