How to solve a nonlinear equation? I'm going to ask my question, explain my problem, show how I thought to solve the problem and where I got stuck, and then re-state my question.
Question: how might I solve this equation:
$$\tag{1} Re_{p}(S)=\frac{S}{1+0.15Re_{p}^{0.687}}$$
for $Re_{p}(S)$?
Explanation of my problem: I am trying to follow a paper, it goes like this:
The settling velocity for a particle falling in a viscous fluid under gravity is
$$\tag{2} v_g=\sqrt{\frac{4}{3} \frac{dg}{C_D} \left(\frac{\rho_p-\rho}{\rho}\right)}$$
where $C_D$ is the drag coefficient, a function of the particle Reynolds number
$$\tag{3} Re_p=\frac{v_g d}{\nu}$$
The correlation for $C_D$ is given as
$$\tag{4} C_D=\frac{24}{Re_p}\left(1+0.15Re_{p}^{0.687}\right)$$
The three equations (2)-(4) can be reduced to a single nonlinear equation with respect to the Reynolds number
$$\tag{5} Re_{p}(S)=\frac{S}{1+0.15Re_{p}^{0.687}}$$
where
$$\tag{6} S=\frac{d^3 g}{18 \nu^2}\left(\frac{\rho_p-\rho}{\rho}\right)$$
The authors of the paper then say the solution of Eqn (5) is shown in the Figure below (as the solid line), which allows one to find the settling velocity
$$\tag{7} v_g=\frac{\nu Re_p(S)}{d}$$
as a function of the settling number $S$.

How I thought to solve the problem and where I got stuck: I thought I might be able to solve for $Re_p(S)$ by rearranging Eqn (1) as so
$$\tag{8} S=Re_p(1+0.15Re_p^{0.687})$$
distributing
$$\tag{9}S=Re_p+0.15Re_p^{1.687}$$
and then procede by performing the "completing the square" method, e.g., dividing through by $0.15$:
$$\tag{10} Re_p^{1.687}+\frac{Re_p}{0.15}=\frac{S}{0.15}$$
and then completing the square by adding $\left(\frac{Re_p}{0.3}\right)^2$ to both sides:
$$\tag{11} Re_p^{1.687}+\frac{Re_p}{0.15}+\left(\frac{Re_p}{0.3}\right)^2=\frac{S}{0.15}+\left(\frac{Re_p}{0.3}\right)^2$$
However, since my first term on the LHS is not to the power of 2 I can't proceed to simplify the LHS.
Back to the question: How did the authors of the paper solve Eqn(1) to produce the solid line in the figure given above?  (Solved for $Re_p(S)$)
(Also, please let me know if my question title is appropriate for the question I am asking)
 A: This is similar to solving equations like $\tan(x) = \log(x)$. There is no way to "isolate" for $x$. In order to solve for the Reynold's numer in your equation, an iterative (i.e. numerical) approach must be taken. I assume you are in engineering, so eventually you will take a Numerical Methods course where you will learn all about solving non-linear equations. 
For now, you can take the following simple iterative approach. Recall your equation
$$ Re_p = \frac{S}{1 + 0.15Re_p^{0.687}}$$
Suppose you want to solve this for $S = 10^0 = 1$. Then,
$$ Re_p = \frac{1}{1 + 0.15Re_p^{0.687}}$$
Make an initial guess for $Re_p$, call it $Re_p^0 = 10$. Use the above equation to iteratively get better guesses for $Re_p$. That is,
$$ Re_p^{i + 1} = \frac{1}{1 + 0.15{Re_p^{i}}^{0.687}}$$
The following table summarizes the results for the first few iterations.
$$\begin{array}{c | c}
i & Re_p^i \\
\hline
0 & 10\\
1 & 0.5772098\\ 
2 & 0.9068355\\ 
3 & 0.8770160\\ 
4 & 0.8794782\\ 
5 & 0.8792734\\ 
6 & 0.8792904\\ 
7 & 0.8792890\\ 
8 & 0.8792891\\ 
9 & 0.8792891\\
\end{array}$$
Of course, you'd want to automate this by writing a simple computer program.
As a final note, I suspect that the authors might have gotten values of $S$ by simply substituting values for $Re_p$ in the above equation, and plotting the result for a range of $Re_p$ that are of interest.
A: $\color{Red}{\textrm{Mathematica code:}}$
R[S_] := Rep /. First@NSolve[S == Rep + 15/100*Rep^(1687/1000), Rep, Reals];

r1[S_] := S; (*  a dummy function *)

LogLogPlot[{R[S], r1[S]}, {S, 10^-1, 10^3}, 
PlotRange -> {Automatic, {10^-1, 10^3/2}}, AxesLabel -> {"Rep", "S"}, 
PlotLabel -> "Settlement of spherical particle in viscous fluid", 
PlotStyle -> {Blue, {Red, Dashed}}, PlotLegends -> {Automatic}]


$\color{Red}{\textrm{Maple code:}}$
restart:
R := proc (S) fsolve(S = Rep+0.15*Rep^1.687, Rep, maxsols = 1) end proc:
plots:-loglogplot(['R(S)', S], S = .1 .. 10^3, view = [0.1 .. 1000, 0.1 .. 500]);

$\color{Red}{\textrm{Maxima code:}}$   
R(S):= find_root (S=Rep + 0.15*Rep^1.678, Rep, 0.10,999)$;
plot2d ([R(S),S], [S, 0.11, 999],logx,logy)$

A: As already said, you need some numerical method to solve the nonlinear equation
$$S=R+k\,R^{\,a}\qquad \text{with}\qquad k=0.15\qquad \text{and}\qquad a=1.687$$ and, as usual, if you want to save iterations, you need a "good" guess for starting.
Playing with inequalities, you know that the solution is somewhere between
$$R_1=\frac S {1+k}\qquad \text{and}\qquad R_2=\left(\frac S {1+k}\right)^{1/a}$$
So, let us take as estimate the geometric mean of these bounds $$R_0=\sqrt{R_1 R_2}=\left(\frac S {1+k}\right)^{\frac{1+a}{2a}}$$
Now, consider that you look for the zero of the function
$$f(R)=R+k\,R^{\,a}-S$$ $$f'(R)=1+a\, k\, R^{\,a-1}$$ Newton iterates will be given by
$$R_{n+1}=R_n-\frac{R_n+k\, R_n^{\,a}-S}{1+a\, k\, R_n^{\,a-1}}$$
For illustration purposes, let us try with $S=10^3$; Newton method will provide the following iterates
$$\left(
\begin{array}{cc}
 n & x_n \\
 0 & 219.186 \\
 1 & 170.109 \\
 2 & 165.981 \\
 3 & 165.950
\end{array}
\right)$$ which is the solution for six significant figures.
For sure, this was for computing a single value. If you want to build the curve, for a second point, use the previous result as estimate. Suppose that we want now $S=1100$; then, the iterates will be
$$\left(
\begin{array}{cc}
 n & x_n \\
 0 & 165.950 \\
 1 & 176.500 \\
 2 & 176.303
\end{array}
\right)$$
For example, using $S=1$, Newton iterates would be
$$\left(
\begin{array}{cc}
 n & x_n \\
 0 & 0.894667 \\
 1 & 0.879285 \\
 2 & 0.879268
\end{array}
\right)$$
Edit
You could do simpler (not to say much better) if you try to inverse the model. Generating $S$ from given $R$ and using least-square method to minimize the sum of squares of relative errors, with a very simplistic model, you could end with
$$R\approx 1.35463\, S^{0.869369}-0.384677\, S$$ which would give pretty good guesses to be used for $R_0$.
For the worked examples, the values generated as estimates are $164.772$, $173.771$ and $0.969949$. Newton method will probably require a couple of iterations for convergence to high accuracy.
