Question: Suppose that $C$ and $B$ are $n \times p$ and $p \times n$, respectively, and that $A + BC$ and $A$ are full rank. Then why is $(I + CA^{-1}B)$ also full rank?

I've tried to manipulate things around, and I wanted to show that 
$(I + CA^{-1}B)x = 0$ implies $x =0$. So far, all I can see is that if $(I + CA^{-1}B)x = 0$, then $(IB + BCA^{-1}B)x = 0$ and thus $(A + BC)(A^{-1}B)x = 0$, so $A^{-1}Bx = 0$ by the non singularity of $A + BC$, and hence $x \in \mathcal{N}(B)$. But that isn't good enough, since $B$ is not full rank.

Answer: 
If you consider using Newton method to find the zero of
$$f(x)=x^3+x^2 (\text{Cb}+\text{Ka})-x (\text{Ca} \,\text{Ka}+\text{Kw})-\text{Ka}
 \, \text{Kw}$$
$$f'(x)=3 x^2+2 x (\text{Cb}+\text{Ka})-(\text{Ca} \text{Ka}+\text{Kw})$$
$$f''(x)=6 x+2 (\text{Cb}+\text{Ka})$$ you must not start at $x=0$ since
$$f(0) \,f''(0)=-2\, \text{Ka} \,\text{Kw} \,(\text{Cb}+\text{Ka})<0$$ (think about Darboux-Fourier theorem). If you use $x\_0=0$, then the first iterate would be
$$x\_1=-\frac{\text{Ka} \, \text{Kw}}{\text{Ca} \,\text{Ka}+\text{Kw}} <0$$

The first derivative cancels at
$$x\_\*=\frac{1}{3} \left(\sqrt{(\text{Cb}+\text{Ka})^2+3 (\text{Ca}\,
 \text{Ka}+\text{Kw})}-(\text{Cb}+\text{Ka})\right) > 0$$
 and the second derivative test shows that this is a minimum.

In order to generate the starting point $x\_0$, perform around $x=x\_\*$ (were $f'(x\_\*)=0$), the Taylor expansion
$$0=f(x\_\*)+\frac 12 f''(x\_\*)\,(x-x\_\*)^2+O((x-x\_\*)^3 ) \implies x\_0=x\_\*+ \sqrt{-2\frac {f(x\_\*)}{f''(x\_\*)}}$$

Using your numbers, Newton iterates would be
$$\left(
\begin{array}{cc}
 n & x\_n \\
 0 & \color{blue}{8.8902}609452502354578 \times 10^{-6}\\
 1 & \color{blue}{8.89021155}71665711819 \times 10^{-6}\\
 2 & \color{blue}{8.8902115568921917511} \times 10^{-6}
\end{array}
\right)$$ This is a very reliable procedure.

**Edit**

I am *almost* sure that we could even skip the Newton procedure using a $[1,2]$ Padé approximant built around $x\_1$. This would give
$$x\_1=x\_0+\frac 12\frac{ f(x\_0) \left(f(x\_0)\, f''(x\_0)-2\, f'(x\_0)^2\right)}{f(x\_0)^2 +f'(x\_0)^3- f(x\_0)\,f'(x\_0)\, f''(x\_0)}$$

For the worked example, this would give $x\_2=\color{blue}{8.890211556892191751}2 \times 10^{-6}$
According to the [computation on Wolfram Alpha](http://www.wolframalpha.com/input/?i=solve+x%5E3+%2B+(0.2%2B1e-4.75)*x%5E2++-(0.1*1e-4.75+%2B+1e-14)*x+%E2%88%921e-4.75*1e-14+%3D+0) linked in the question, the task is to find the roots of the cubic equation 

$x^{3} + (0.2 + 1 \times 10^{-4.75})x^{2}-(0.1 \times 10^{-4.75} + 1 \times 10^{-14})x - (1 \times 10^{-4.75} \times 1 \times 10^{-14})$

If this is representative of the cubic equations that need to be solved, there seems to be no problem in tackling this with IEEE-754 double precision computation and a standard solver, such as Kahan's `QBC` solver:

William Kahan, "To Solve a Real Cubic Equation". Lecture notes, UC Berkeley, Nov. 1986 ([online](https://cs.berkeley.edu/~wkahan/Math128/Cubic.pdf))

I am using Kahan's code pretty much verbatim, but have enhanced it in various places by the use of the FMA (fused multiply-add) operation which is readily accessible from C++ as the standard math function `fma()`. I specified the coefficients as follows:

```
double a3 =  1.0;
double a2 =  0.2 + exp10 (-4.75);
double a1 = -(0.1 * exp10 (-4.75) + exp10 (-14.0));
double a0 = -(exp10 (-4.75) * exp10 (-14.0));

```

The results returned by the `QBC` solver are:

```
-2.0002667300555729e-001
-9.9999998312876059e-014 + i *  0.0000000000000000e+000
 8.8902115568921925e-006 + i *  0.0000000000000000e+000

```

These match the results computed by Wolfram Alpha: $-0.200027, -1.\times10^{-13}, 8.89021\times10^{-6}$