The two circles $C_r$ and $C_s$ are mapped to the parrallel lines $\left\{ \operatorname{Re} w = \frac{1}{2r}\right\}$ and $\left\{ \operatorname{Re} w = \frac{1}{2s}\right\}$ by the inversion. Since the inversion maps circles to circles (where a straight line is counted as a circle of infinite radius), it maps each of the $C_i$ to a circle touching both of these parallel lines. Also, the image of $C_1$ touches the real axis (since $C_1$ does, and the inversion maps the real axis $\cup \{\infty\}$ to itself), and further, each circle $C_k$ for $k > 1$ touches the circles $C_{k-1}$ and $C_{k+1}$, hence so do their images.
The distance between the two parallel lines is $\frac{1}{2s} - \frac{1}{2r}$, so the images of the $C_k$ have the common radius
$$R = \frac{1}{2}\left(\frac{1}{2s} - \frac{1}{2r}\right),$$
and their centres all lie on the line $$\left\{ \operatorname{Re} w = M\right\};\qquad M := \frac{1}{2}\left(\frac{1}{2s} + \frac{1}{2r}\right).$$
It follows that the centre of the image $C_k'$ of $C_k$ is
$$b_k = M - (2k-1)\cdot iR,$$
and hence the defining equation of $C_k'$ is
$$\lvert w-b_k\rvert^2 = R^2,$$
or
$$w\overline{w} - b_k\overline{w} - \overline{b_k}w + (\lvert b_k\rvert^2 - R^2) = 0.$$
Applying the inversion, we see that $C_k$ has the defining equation
$$1 - b_k z - \overline{b_kz} + (\lvert b_k\rvert^2-R^2)z\overline{z} = 0.$$
Since $\lvert b_k\rvert^2 = M^2 + (2k-1)^2R^2 > R^2$, we can divide and obtain the equivalent equation
$$z\overline{z} - \frac{\overline{b_k}}{\lvert b_k\rvert^2-R^2}\overline{z} - \frac{b_k}{\lvert b_k\rvert^2 - R^2} z + \frac{1}{\lvert b_k\rvert^2-R^2} = 0,$$
or
$$\left\lvert z - \frac{\overline{b_k}}{\lvert b_k\rvert^2-R^2} \right\rvert^2 = \left(\frac{R}{\lvert b_k\rvert^2-R^2}\right)^2.$$
The desired relation between the radii should not be difficult to obtain from that.
Most programs that support Matlab-style or NumPy-style array operations and plotting can do that. Examples are Octave, FreeMat, or Scilab.
The essence (in Octave syntax), demonstrated for the $\tanh$ function, is:
n=30;
[X,Y]=ndgrid([-pi/2:pi/n:pi/2],[-pi/2:pi/n:pi/2]);
Z=X+1i*Y;
W=tanh(Z);
clf(); axis([-4,4,-3,3],"equal"); hold on;
plot(W);
plot(W');
But then the segments between the grid points will be straight. You may want a little more refinement. Run the following as a script file or enter its contents in the Octave UI:
#! /usr/bin/env octave -qf
# Set plot ranges and aspect ratio
clf(); axis([-4,4,-3,3],"equal"); hold on;
# Set plot title
title("Map of tanh(z) resp. isolines of artanh(z)");
# Set number of subintervals of the domain used (in each direction)
n=360;
# Plot only each m-th curve. This makes hi-res arcs appear curved.
m=12;
# The square domain, subdivided
[X,Y]=ndgrid([-pi/2:pi/n:pi/2],[-pi/2:pi/n:pi/2]); Z=X+1i*Y;
# The image
W=tanh(Z);
# Plot image of every m-th line for which real(Z)=const.
# Hint: Specify ";;" to suppress the legend for octave < 2.9.13
plot(real(W(1:m:n+1,:))',imag(W(1:m:n+1,:))',";;");
# Plot image of every m-th line for which imag(Z)=const.
plot(real(W(:,1:m:n+1)),imag(W(:,1:m:n+1)),";;");
print("-dpng", "-mono", "-solid", "-S480,360", "tanh-map.png");
Best Answer
https://github.com/search?l=C%2B%2B&q=complex+functions&ref=cmdform&type=Repositories
http://www.codeproject.com/Articles/80641/Visualizing-Complex-Functions
http://webs.um.es/gvb/AC/FC.zip
Maybe you find what you're looking in the first two links, the third is in Spanish and need Ms-Dos to run, I remember that was possible choosing strips.
But without doubt the best method is programming directly with help of some graphics library and feel the freedom of making what you need at any moment. I use to work with C and OpenGL and works really fast. You can insert parameters, modify them in real time and make what you want. That's pretty fun.