Special conformal killing fields - solving for integral curves. For each $b\in\mathbb R^d$, let a vector field $X_b:\mathbb R^d\to\mathbb R^d$ be defined as follows:
\begin{align}
  X_b(x) = 2(b\cdot x)x - x^2 b,
\end{align}
where $x^2 = x\cdot x$.  This is the conformal killing field that generates special conformal transformations.

How does one solve for the integral curves of this monster?

Apparently, the (local) answer (which I've only ever seen magically pulled out of a hat) is as follows:
\begin{align}
  x(t) = \frac{x_0 - x_0^2(tb)}{1-2x_0\cdot(tb) + x_0^2(tb)^2} \tag{$\star$}
\end{align}
How can one at least somewhat systematically arrive at this?  The system of ODEs one needs to solve is non-linear and coupled; I'm rather at a loss as to how to even begin attacking it.
Note: eq. $(\star)$ written down and discussed briefly in 
A Mathematical Introduction to Conformal Field theory by Schottenloher, 2nd ed., p. 17.
It is also all over the place in the theoretical physics literature where one sees statements to the effect of "integrating the infinitesimal version of the special conformal transformation gives..." and then $(\star)$ is written down.
Progress! I think I've almost figured it out using a tricky change of variables.  We want to solve
\begin{align}
  \dot x = 2(b\cdot x)x - x^2 b.
\end{align}
Make a change of variables
\begin{align}
  y = \frac{x}{x^2}
\end{align}
and after some algebra, one shows that $y$ satisfies
\begin{align}
  \dot y = -b
\end{align}
which has as a solution
\begin{align}
  y = y_0 - tb,
\end{align}
so now we simply need to solve the algebraic equation
\begin{align}
  \frac{x}{x^2} = \frac{x_0}{x_0^2} - tb.
\end{align}
The solution quoted above certainly solves this, but how does one solve this equation "from scratch?"
Solution As pointed out by, @HunsLundmark, the algebraic equation
\begin{align}
  \frac{x}{x^2} = A
\end{align}
can be readily solved to give $x = A/A^2$, so setting $A = x_0/x_0^2-tb$ gives the desired result.
 A: 1. Introduce a new independent variable $s=bt$ to get rid of $b$, and
w.l.o.g., assume that the vector $b$ points along $e_1$, $b=be_1$:
$$ \frac{dx}{ds} = 2x_1x-x^2e_1, $$
$$ \dot x_1 = x_1^2-x_2^2-\cdots-x_n^2, \qquad \dot x_i = 2x_1x_i, $$
where, since $x_1$ is special, the index $i$ ranges over $2,\ldots,n$.
2. Introduce the change of variables $(z, w) = \left(x_1, (\sum_i
  x_i^2)^{1/2}\right)$, so that we need to solve
$$ \dot z = z^2-w, \qquad \dot w = 4zw. $$
Differentiate $\dot z$ to get
$$ \ddot z = 6z\dot z-4z^3. $$
3. Now here I cheated and used Mathematica: it gives the general solution as
$$ z(s) = -\frac{1}{2}\left(\frac{1}{s+\gamma} +
  \frac{1}{s+\bar\gamma}\right), $$
where $\gamma = \alpha + i\beta$ is a complex constant.
3a. Another way is
to look for poles of $z$ by substituting $z(s)=a(s+\gamma)^b$,
expanding in Laurent series at $s=-\gamma$ and observing that the
resulting leading order equation
$$ a b(b-1)(s+\gamma)^{b-2} =
6a^2b(s+\gamma)^{2b-1}-4a^3(s+\gamma)^{3b} $$
has solutions $b=-1$ and $a=0, -\frac12, -1$. Also note that the
Laurent series at $s=\infty$ also gives $z\sim -(\text{$1$ or $\frac12$})s^{-1}$. A complex-differentiable function with only regular singularities on $\mathbb{C}\cup\{\infty\}$ must be rational, so only $a=-1$ with one pole and $a=-\frac12$ with two poles make sense. A rational function with two order-1 poles and finite at $\infty$ must be the sum of only the two pole terms, as above. Since $z$ should be real, the second pole is $\bar\gamma$. Not completely rigourous, but good enough.
4. Integrating $w$ is easy:
$$ (\log w)' = 4z, \qquad w(s) = w(0)\left(
\frac{\alpha^2+\beta^2}{(s+\alpha)^2+\beta^2} \right)^2. $$
Substituting $z, w$ into $\dot z$ gives the relationship between
$\alpha$ and $\beta$:
$$ \frac{\alpha^2-\beta^2}{(\alpha^2+\beta^2)^2} = z_0^2-w_0, $$
so therefore
$$ \gamma = \frac{-1}{z_0+i w_0}, \qquad \Re\gamma =
\frac{-z_0}{z_0^2+w_0^2}. $$
Integrating $x_i$ is as easy as $w$:
$$ x_i(s) = x_i(0)
\frac{\gamma\bar\gamma}{(s+\gamma)(s+\bar\gamma)}. $$
5. Finally, writing it out in full, letting $\rho_0^2 =
x_1(0)^2+\cdots+x_n(0)^2$:
$$\begin{aligned}
 \vec x &= \frac{1}{(s+\gamma)(s+\bar\gamma)}
\left(-(s+\Re\gamma)e_1+x_i(0)\gamma\bar\gamma \right)
\\&= \frac{-(s\rho_0^2-z_0)e_1 +
  x_i(0)}{\rho_0^2(s+\gamma)(s+\bar\gamma)}
\\&= \frac{\vec x(0) - s\rho_0^2e_1}{1-2z_0s+\rho_0^2s^2}.
\end{aligned}$$
This is clearly identical to your formula, once $t=s/b$ and $b=be_1$
are substituted back in.
