An affine transformation $az+b$ may be used to turn one circle into the unit circle (specifically, $(z-a)/r$ if the circle has center $a$ and radius $r$). If the other circle is not contained within it, we may apply $1/x$ to (or $-1/x$ if we want to use $\mathrm{PSL}_2\mathbb{R}$), which fixes the unit circle while swapping its interior and exterior.
The transformations that preserve the real line come from $\mathrm{PSL}_2\mathbb{R}$ and the ones that preserve the unit circle are $\mathrm{PSU}(1,1)$. The intersection of these two is $\mathrm{PSO}(1,1)$ whose elements may be represented by
$$ \begin{bmatrix} \cosh\tau & \sinh\tau \\ \sinh\tau & \cosh\tau \end{bmatrix}. $$
Alternatively, one may say $Mz=\tfrac{az+b}{cz+d}$ and solve $M(1)=1$, $M(-1)=-1$ to get $Mz=\frac{az+b}{bz+a}$, without loss of generality $a^2-b^2=1$ whose solutions ($a>0$) are parametrized $(a,b)=(\cosh\tau,\sinh\tau)$.
The circle in the interior of the unit circle intersects at points $u,v\in(-1,1)$. We may parametrize this interval with $\tanh$, say $u=\tanh(\alpha)$, $v=\tanh(\beta)$ and check that
$$ \begin{bmatrix} \cosh\tau & \sinh\tau \\ \sinh\tau & \cosh\tau \end{bmatrix} \tanh(\sigma)=\tanh(\tau+\sigma) $$
using $\tanh$'s sum formula. Therefore, we may choose $\tau=-\frac{1}{2}(\alpha+\beta)$ and apply this transformation to turn $u,v$ into values $x,-x$ within $(-1,1)$, necessarily intersection points of a concentric circle.
Note the transformation may also be written as
$$ \begin{bmatrix} 1 & \tanh\tau \\ \tanh\tau & 1 \end{bmatrix} $$
and use various hyperbolic trig identities to obtain
$$ \tanh\tau=\pm\sqrt{\frac{1+uv-\sqrt{(1-u^2)(1-v^2)}}{1+uv+\sqrt{(1-u^2)(1-v^2)}}} $$
In the case of $C(0,1)$ and $C(4,2)$, we can apply $-1/x$ so the latter becomes $C(-\tfrac{1}{3},-\tfrac{1}{6})$ with corresponding endpoints $u=-\tfrac{1}{2}$, $v=-\tfrac{1}{6}$. Plugging these into the above, rescaling and using Wolfram|Alpha to denest,
$$ \begin{bmatrix} \sqrt{21}+\sqrt{5} & \sqrt{21}-\sqrt{5} \\ \sqrt{21}-\sqrt{5} & \sqrt{21}+\sqrt{5} \end{bmatrix} $$
preserves $C(0,1)$ and turns $C(-\tfrac{1}{3},-\tfrac{1}{6})$ into $C(0,\tfrac{11-\sqrt{105}}{4})$.
Best Answer
The real axis is a line of symmetry. So I expect that, to centre the new circles on the origin, we need $f(1/2)=-f(0)$ and $f(-1)=-f(1)$. So $$(a+2b)/(c+2d)=-b/d\\(a+b)/(c+d)=-(a-b)/(c-d)$$
Pick almost any values for $a$ and $c$, and solve for $b$ and $d$
for example, $a=c=1$, then $$(1+2b)d=-b(1+2d),(1+b)(1-d)=-(1-b)(1+d)\\ b+d+4bd=0,2-2bd=0\\ b+d=-4,bd=1\\ b^2+4b+1=0$$