Main Question
I am following the book Nonlinear Dynamics: A concise introduction interlaced with code, and in chapter 2, there is a discussion about plotting nullclines, with the example of the Fitzhugh-Nagumo model
\begin{equation} \begin{aligned} \dot{u} &= au(u-b)(1-u) - w + I\\ \dot{w} &= \varepsilon(u - w) \end{aligned} \end{equation}
where a = 3, b = 0.2, I = 0, and $\varepsilon = 0.01$. I understand how for plotting nullclines you are interested in where the individual state variables' rate of change becomes zero, but in the codebase used to produce a figure 2.4 (see the official github), they define the nullclines as below:
\begin{equation} \begin{aligned} \dot{u} = 0 &\Leftrightarrow n_u(u) = au(u-b)(1-u) \\ \dot{w} = 0 &\Leftrightarrow n_w(u) = u \end{aligned} \end{equation}
I see how
\begin{equation} \begin{aligned} \dot{w} &= 0 \\ 0 &= \varepsilon(u - w) \\ w &= u \\ n_w(u) &= u \end{aligned} \end{equation}
but how does $n_u(u)$ have only a single solution? For example and using the aforementioned parameter values,
\begin{equation} \begin{aligned} \dot{u} &= 0 \\ 0 &= au(u - b)(1 - u) - w \\ u &=\ ... && \text{Solve[au(u -b)(1-u) - w == 0, u]} \end{aligned} \end{equation}
where Solve[au(u -b)(1-u) - w == 0, u] results in multiple solutions (after checking Mathematica's outputs), and such solutions do not match the one provided by the official github. It seems that what the authors of the book do (and this is also apparent on the wiki page for Fitzhugh-Nagumo model), is to just isolate $w$, and then $w$ just becomes a function of $u$.
Context
I am trying to plot nullclines in figure S1 of "Quantitative modeling of the terminal differentiation of B cells and mechanisms of lymphomagenesis", and the ODEs for the desired plot are defined below but would also end up with multiple roots for $\dot{b}$.
\begin{equation} \begin{aligned} \dot{p} &= \mu_p + \sigma_p \frac{k_b^2}{k_b^2 + b^2} - \lambda_p p \\ \dot{b} &= \mu_b + \sigma_b \frac{k_p^2}{k_p^2 + p^2}\frac{k_b^2}{k_b^2 + b^2} - (\lambda_b + bcr_0\frac{k_b^2}{k_b^2 + b^2})b \end{aligned} \end{equation}
Note that $\mu_p$, $\mu_b$, $\sigma_p$, $\sigma_b$, $k_b$, $k_p$, $\lambda_b$, $\lambda_p$, and $bcr_0$ are all parameters of the model. I can also easily see that
\begin{equation} \begin{aligned} \dot{p} &= 0 \\ 0 &= \mu_p + \sigma_p \frac{k_b^2}{k_b^2 + b^2} - \lambda_p p \\ \lambda_p p &= \mu_p + \sigma_p \frac{k_b^2}{k_b^2 + b^2} \\ p &= \frac{\mu_p + \sigma_p \frac{k_b^2}{k_b^2 + b^2}}{\lambda_p} \\ n_p(b) &= \frac{\mu_p + \sigma_p \frac{k_b^2}{k_b^2 + b^2}}{\lambda_p} \end{aligned} \end{equation}
but for $\dot{b} = 0$, I am not so certain since multiple roots will result from Solve[$\dot{b}$ == 0, b], and if I were to try to instead isolate $p$, I end up with the below expression (which when plotting the nullclines for $\dot{p}$ and $\dot{b}$ does match figure S2). Here is my solution for the nullcline of $\dot{b}$ along with my plot first and plot from the paper:
\begin{equation} \begin{aligned} \frac{db}{dt} &= \mu_b + \sigma_b \frac{k_p^2}{k_p^2 + p^2}\frac{k_b^2}{k_b^2 + b^2} - (\lambda_b + BCR)b \\ \frac{db}{dt} &= 0 \\ \mu_b + \sigma_b \frac{k_p^2}{k_p^2 + p^2}\frac{k_b^2}{k_b^2 + b^2} - (\lambda_b + BCR)b &= 0 \\ \mu_b + \sigma_b \frac{k_p^2}{k_p^2 + p^2}\frac{k_b^2}{k_b^2 + b^2} &= (\lambda_b + BCR)b \\ \sigma_b \frac{k_p^2}{k_p^2 + p^2}\frac{k_b^2}{k_b^2 + b^2} &= (\lambda_b + BCR)b - \mu_b \\ \frac{k_p^2}{k_p^2 + p^2}(\frac{\sigma_b k_b^2}{k_b^2 + b^2}) &= \lambda_b b + (BCR)b - \mu_b \\ \frac{k_p^2}{k_p^2 + p^2} &= [\lambda_b b + (BCR)b - \mu_b](\frac{k_b^2 + b^2}{\sigma_b k_b^2}) \\ \frac{k_p^2}{[\lambda_b b + (BCR)b - \mu_b](\frac{k_b^2 + b^2}{\sigma_b k_b^2})} &= k_p^2 + p^2 \\ \sqrt{\frac{k_p^2}{[\lambda_b b + (BCR)b - \mu_b](\frac{k_b^2 + b^2}{\sigma_b k_b^2})} - k_p^2} &= p = n_b(b) \end{aligned} \end{equation}
