In order to determine the finite SCT from its infinitesimal version, we need to solve for the integral curves of the special conformal killing field $X$ defined by
\begin{align}
X(x) = 2(b\cdot x) x - x^2 b.
\end{align}
I explain why this is equivalent to "integrating" the infinitesimal transformation below. This means we need to solve the differential equation $X(x(t)) = \dot x(t)$ for the function $x$. Explicitly, this differential equation is
\begin{align}
\dot x = 2(b\cdot x) x - x^2 b.
\end{align}
This can be done with a trick, namely a certain change of variables. Define
\begin{align}
y = \frac{x}{x^2},
\end{align}
then the resulting differential equation satisfied by $y$ becomes simple;
\begin{align}
\dot y = -b.
\end{align}
I urge you to perform the algebra yourself to confirm this. It's kind of magic that it works if you ask me, and the change of variables is precisely an inversion, so I think there's something deeper going on here, but I'm not sure what it is. The solution to this equation is simply $y = y_0 -tb$, so we find that the original function $x$ satisifes
\begin{align}
\frac{x}{x^2} = \frac{x_0}{x_0^2} - tb.
\end{align}
In other words, we've turned a monstrous nonlinear system of ODEs into a simple algebraic equation. In fact, one can show that the algebraic eqution $x/x^2 = A$ has the solution $x = A/A^2$, from which it follows that the solution to our original problem is
\begin{align}
x(t) = \frac{x_0 - x_0^2(tb)}{1-2x_0\cdot(tb) + x_0^2(tb)^2},
\end{align}
as desired, since this is precisely the form of the "finite" SCT. Note that these only are local integral curves; the solution hits a singularity when $t$ is such that the denominator vanishes.
Why solve for integral curves?
If you're wondering what your original question has anything to do with solving the integral curves of the special conformal vector field I wrote down, then read on.
It helps to start with the concept of a flow.
Transforming points via flows.
Let a point $x\in\mathbb R^d$ be given, then for each $b\in\mathbb R^d$, we assume, at least in a neighborhood of that point, that there is a $\epsilon$-parameter family of transformations $\Phi_b(\epsilon):\mathbb R^d \to \mathbb R^d$ with $\epsilon\in [0,\bar\epsilon)$ such that $\Phi_b(\epsilon)(x)$ tells you what an SCT corresponding to the vector $b$ does to the point $x$ is you ``flow" in $\epsilon$. At $\epsilon = 0$, this flow just maps the point to itself;
\begin{align}
\Phi_b(0)(x) = x,
\end{align}
namely it starts at the identity. For $\epsilon >0$, the flow translates the point along a curve in $\mathbb R^d$. If you change $b$, then this corresponds to moving way from $x$ in a different initial direction under the flow.
Infinitesimal generator of a flow.
When we talk of an infinitesimal generator of such a transformation, we are talking about the term that generates the linear approximation for the flow in the parameter $\epsilon$. In other words, we expand
\begin{align}
\Phi_b(\epsilon)(x) = x + \epsilon G_b(x) + O(\epsilon^2),
\end{align}
and the vector field $G_b$ is called the infinitesimal generator of the flow. What you have pointed out in your question is that we know this infinitesimal generator;
\begin{align}
G_b(x) = 2(b\cdot x)x - x^2 b,
\end{align}
and we now want to reconstruct the whole flow simply by knowing this information corresponding to its linear approximation at every point. This is equivalent to solving some first order ordinary differential equations, which is why people often say we want to "integrate" the infinitesimal transformation to determine the finite one; integration is a perhaps somewhat archaic way of solving the corresponding differential equation.
Finding the whole flow given its generator.
Ok, so what differential equation do we solve? Well, note that the vector field $G_b$ is tangent to the curves generated by the flow by its very construction (we took a derivative with respect to $\epsilon$ with is the "velocity" of the curve generated by the flow), so the differential equation we want to solve is
\begin{align}
\dot x(\epsilon) = G_b(x(\epsilon)),
\end{align}
and we want to solve for $x(\epsilon)$. The solutions to this differential equation are referred to as integral curves of the vector field $G_b$.
Acknowledgements.
I didn't figure out the first part of this answer completely on my own. The idea for making the substitution $y=x/x^2$, which is really the crux of everything, came from here http://www.physicsforums.com/showthread.php?t=518316, namely from user Bill_K.
The idea for how to solve the algebraic equation $x/x^2 = A$ came from math.SE user @HansLundmark after I posted essentially your question in mathy language on math.SE here.
I should another math.SE user @Kirill solved for the integral curves in a totally different way in his answer to the question I posted.
Addendum.
How does one get $\dot y = -b$ from the change of variable $y=x/x^2$ as claimed in the first section? Well, let's calculate:
\begin{align}
\dot y
&= \frac{x^2\dot x - x(2x\cdot \dot x)}{(x^2)^2} \\
&= \frac{x^2(2(b\cdot x) x - x^2 b) - x(2x\cdot (2(b\cdot x) x - x^2 b))}{(x^2)^2} \\
&= \frac{2x^2(b\cdot x) x - (x^2)^2b - 4x^2(b\cdot x) x + 2x^2(b\cdot x)x}{(x^2)^2} \\
&= -\frac{(x^2 )^2b}{(x^2)^2} \\
&= -b
\end{align}
Magic!
The fields transform under finite conformal transformations as$^1$
$$
\Phi^a(x') \mapsto {\Phi^a}'(x) = \Omega(x')^\Delta\,D(R(x'))^{\phantom{b}a}_b \,\Phi^b(x')\,. \tag{1}\label{main}
$$
as given in equation $(55)$ of $[1]$. Let's break it down:
- $\Delta$ is the conformal dimension of $\Phi$.
- $\Omega$ is the conformal factor of the transformation.
- $D$ is the spin representation of $\Phi$.
- $R$ is the rotation Jacobian of the transformation
So let's compute these things. The spin and the conformal dimension $\Delta$ are given. The first thing we have to look at is the Jacobian.
$$
\frac{\partial x^{\prime \mu}}{\partial x^\nu} = \Omega(x') R^\mu_{\phantom{\mu}\nu}(x')\,.
$$
This implicitly defines both $\Omega$ and $R$ and it is not ambiguous because we require $R \in \mathrm{SO}(d)$, namely
$$
R^{\mu}_{\phantom{\mu}\nu} \,\eta^{\nu\rho}\,R^{\lambda}_{\phantom{\lambda}\rho} \,\eta_{\lambda\kappa}= R^{\mu}_{\phantom{\mu}\kappa}\,.
$$
You can immediately see that for the Poincaré subgroup of the conformal group $\Omega(x')= 1$, whereas for dilatations $\Omega(x') = \lambda$ and for special conformal transformations
$$
\Omega(x') = \frac{1}{1+2(b\cdot x') + b^2 {x'}^2}\,.\tag{2}\label{omega}
$$
This can be proven with a bit of algebra. Using your definition of SCT
$$
{x'}^\mu = \frac{x^\mu - b^\mu x^2}{1+-2(b\cdot x) + b^2 {x}^2}\,,
$$
one can check
$$
\frac{\partial x^{\prime \mu}}{\partial x^\rho} \eta_{\mu\nu} \frac{\partial x^{\prime \mu}}{\partial x^\lambda} = \frac{\eta_{\rho\lambda}}{(1-2(b\cdot x) + b^2 x^2)^2}\,.
$$
That means that the Jacobian is an orthogonal matrix up to a factor, which is the square root of whatever multiplies $\eta_{\rho\lambda}$. Then we have to re-express that as a function of $x'$. After some algebra again one finds that it suffices to change the sign to the term linear in $b$.
Finally, how does one compute $R$? Well, it's just the Jacobian divided by $\Omega$. In the case of special conformal transformations one has (there might be mistakes, redo it for safety)
$$
R^{\mu}_{\phantom{\mu}\nu} = \delta^\mu_\nu + \frac{2 b_\nu x^\mu - 2 b^\mu (b_\nu x^2+ x_\nu - 2 (b\cdot x) x_\nu) -2b^2 x^\mu x_\nu }{1-2b\cdot x +b^2 x^2}\,,
$$
which, as before, needs to be expressed in terms of $x'$.
If you are interested in $\Phi$ scalar then $D(R) = 1$ and you can just plug \eqref{omega} into \eqref{main} to obtain the transformation. If you want to consider also spinning $\Phi$ then it's not much harder.
For spin $\ell=1$ the $D$ is just the identity, namely
$$
D(R)^{\phantom{\nu}\mu}_\nu = R^{\phantom{\nu}\mu}_\nu\,.
$$
For higher spins one just has to take the product
$$
D(R)^{\phantom{\nu_1\cdots \nu_\ell}\mu_1\cdots \mu_\ell}_{\nu_1\cdots \nu_\ell} = R^{\phantom{\nu_1}\mu_1}_{\nu_1}\cdots R^{\phantom{\nu_\ell}\mu_\ell}_{\nu_\ell}\,.
$$
Again, by plugging these definitions in \eqref{main} you obtain the desired result.
$\;[1]\;\;$TASI Lectures on the Conformal Bootstrap,
David Simmons-Duffin, 1602.07982
$\;{}^1\;\;$The way the transformations are written in the lecture notes linked above differs a bit from Di Francesco Mathieu Sénéchal. The difference is that Di Francesco et al. make an active transformation $x \to x'$ with
$$
\Phi(x) \mapsto \Phi'(x') = \mathcal{F}(\Phi(x))\,,
$$
while David Simmons Duffin makes essentially the inverse transformation $x' \to x$
$$
\Phi(x') \mapsto \Phi'(x) = \mathcal{F}^{-1}(\Phi(x'))\,.
$$
That is why in the above discussion the indices of $R^\mu_{\phantom{\mu}\nu}$ get swapped when passed inside $D$ as $D(R) = R^{\phantom{\nu}\mu}_{\nu}$. And that's also why we get a factor $\lambda^\Delta$ rather than $\lambda^{-\Delta}$ as Di Francesco et al. This is all consistent as long as it is clear what one is doing.
Best Answer
@Trimok solved the problem most elegantly in his comment to the question cited, and since you are troubled by @Josh's simplifying changes of variables, $b^\mu\equiv \hat{b}~ t$, I'll avoid them to merely integrate the variation $$\delta x^\mu =2(b\cdot x)x^\mu-x\cdot x~ b^\mu$$ directly. It immediately implies $$\delta \left (\frac{ x^\mu}{x^2}\right) = \frac{ \delta x^\mu}{x^2} -2 x^\mu \frac{ x\cdot \delta x }{x^4} = -b^\mu ~.$$
That is, the vector $x^\mu/x^2$ evolves by shifting along $-\hat{b}^\mu$, linearly in the magnitude $|b|$ of $b^\mu$, so with constant unit speed in this "pseudotime" $|b|$. Integrating this simplest of advections for finite pseudotime, we immediately get $$ \frac{ x'^\mu}{x'^2}= \frac{ x^\mu}{x^2} -b^\mu . $$ Square both sides, to get the normalization, $$ \frac{ 1}{x'^2}= \frac{1}{x^2} +b^2 -2\frac{ x\cdot b}{x^2}= \frac{1-2x\cdot b+b^2 x^2 }{x^2}, $$ which divides the above vector equation to yield your conventional form for it, ∴ $${x^\mu}'=\frac{x^\mu-b^\mu x^2}{1-2x\cdot b+b^2 x^2}~.$$