I am trying to prove this function has a unique positive root in $\lambda$ for any $\beta>0$ and $\eta>0$ (as background, this is the first derivative of a function I try to prove is unimodal). Graphically, I cannot find any counterexample. I can prove the function has a root since $f(0)>0$ and $\lim_{\lambda\rightarrow \infty} f\left(\lambda\right)<0$; but would like to go one step further and show that this root is unique. Any help is much appreciated!
Thanks!
$$ f\left(\lambda \mid \beta, \eta \right)= \frac{1}{\beta +\eta-e^{-\frac{\lambda}{\beta }} (\beta +\lambda) }-\frac{1}{\eta + e^{-\frac{\lambda}{\beta }} (\beta +\lambda)} -\frac{\lambda^2 e^{-\frac{\lambda}{\beta }} \left(\frac{1}{\left(\eta + e^{-\frac{\lambda}{\beta }} (\beta +\lambda) \right)^2}+\frac{1}{\left(\beta +\eta-e^{-\frac{\lambda}{\beta }} (\beta +\lambda) \right)^2}\right)}{\beta }$$
P.S. Strictly speaking, I only need to show uniqueness on $\lambda\in\left[0, -\beta\left(1+W_{-1}\left(-\frac{1}{2e}\right)\right)\right]$, where $W_{-1}$ is the secondary branch of the Lambert function. But my experimenting with the function seems to suggest the root is unique for any positive $\lambda$.
Not a complete answer, but I think it's close with some additional inequality-chasing on your part.
Denote $a(\lambda, \beta, \eta) = \eta + e^{-\lambda/\beta}(\beta+\lambda)$ and $b(\lambda, \beta, \eta) = \beta + \eta - e^{-\lambda/\beta}(\beta+\lambda)$. Since $a(\cdot), b(\cdot) > 0$ for all $\lambda, \beta, \eta > 0$, $f(\lambda|\beta,\eta)$ has a unique root in $\lambda$ if and only if \begin{align*} g(\lambda, \beta, \eta) \overset{\text{def}}{=} a(\lambda, \beta, \eta)b(\lambda, \beta, \eta)f(\lambda|\beta, \eta) = a(\lambda, \beta, \eta) - b(\lambda, \beta, \eta) - \frac{\lambda^2}{\beta}e^{-\lambda/\beta}\left(\frac{a}{b} + \frac{b}{a}\right) \end{align*} has a unique root in $\lambda$. Henceforth, I will sometimes denote functions as $a(\lambda), b(\lambda), g(\lambda)$, or sometimes without arguments as $a, b, g$ if there is no confusion.
First, note that for any $k \in \mathbb{R}$ \begin{align*} g(k\lambda, k\beta, k\eta) = k g(\lambda, \beta, \eta) \end{align*} Showing that $g$ is 1st-order homogeneous. Without loss of generality, we can set $\eta$ to any constant to only deal with two arguments $(\lambda, \beta)$. (In my search, I set $\eta = 1$, but will keep the derivations below as general as possible.)
Secondly, to prove that $g$ has a unique root in $\lambda$, it suffices to prove that for all $\beta > 0$ there exists an $\xi > 0$ such that \begin{align*} g'(\lambda, \beta) < 0, \lambda \in (0, \xi) \quad \text{and} \quad g(\lambda, \beta) < 0, \lambda \in [\xi, \infty) \end{align*} Since this means $g$ is decreasing until it hits a negative value, and then stays negative henceforth.
With this in mind, note that \begin{align*} g(\lambda) = - (\lambda - r_1(\lambda))(\lambda - r_2(\lambda)) \end{align*} where $r_1, r_2$ are the smaller and larger values, respectively, of \begin{align*} \frac{2e^{-\lambda/\beta} \pm \sqrt{4e^{-2\lambda/\beta}-4\frac{e^{-\lambda/\beta}}{\beta}(\frac{a}{b}+ \frac{b}{a})(1-2\exp(-\lambda/\beta))\beta}}{2\frac{e^{-\lambda/\beta}}{\beta}(\frac{a}{b} + \frac{b}{a})} \end{align*} If $r_1, r_2$ were not covarying with $\lambda$, then they could be viewed as roots of the quadratic, in which case $g<0$ always if the discriminant (quantity under the square root) is negative. Using the inequality $\frac{a}{b} + \frac{b}{a} \ge 2$, one can show that the discriminant is negative whenever $\lambda > \log(\frac{5}{2})\beta \approx 0.916\beta$. Perhaps we wish this could be our value of $\xi$, but it's not strong enough through my attempts, but at least we can now focus on regions $\lambda \le 0.916\beta$.
For $\lambda \le 0.916\beta$, you can show that $\lambda > r_1(\lambda)$ always, so a necessary and sufficient condition for $g(\lambda) < 0$ is for \begin{align*} \lambda &> r_2(\lambda) \\ &=\frac{2e^{-\lambda/\beta} \pm \sqrt{4e^{-2\lambda/\beta}-4\frac{e^{-\lambda/\beta}}{\beta}(\frac{a}{b}+ \frac{b}{a})(1-2\exp(-\lambda/\beta))\beta}}{2\frac{e^{-\lambda/\beta}}{\beta}(\frac{a}{b} + \frac{b}{a})} \\ &= \beta\left(\left(\frac{a}{b} + \frac{b}{a}\right)^{-1} + \sqrt{\left(\frac{a}{b} + \frac{b}{a}\right)^{-2} + (2e^{-\lambda/\beta}-1)e^{-\lambda/\beta}\left(\frac{a}{b} + \frac{b}{a}\right)^{-1}}\right) \end{align*} and the derivative of $g$ is \begin{align*} g'(\lambda) = 2e^{-\lambda/\beta}\frac{\lambda}{\beta}\left(1 - \frac{\lambda}{\beta}\right) - (\beta + 2\eta)^2 \frac{\lambda}{\beta}e^{-\lambda/\beta}\left(\left(2-\frac{\lambda}{\beta}\right)h(\lambda) - \lambda h'(\lambda)\right)\frac{1}{h^2} \end{align*} where \begin{align*} h(\lambda) &= \eta(\beta+\eta) + \beta e^{-\lambda/\beta}(\beta+\lambda) - e^{-2\lambda/\beta}(\beta+\lambda)^2 \\ h'(\lambda) &= \frac{\lambda}{\beta}e^{-2\lambda/\beta}\left(2\lambda - \beta\left(e^{\lambda/\beta}-2\right)\right) \end{align*} There is a lot of flexibility in finding this function $\xi$; in my simulations, as well as the plot you show, $f$ (and also $g$) becomes negative before it is increasing again, so you just need to find a function $r^*(\beta)$ between those two two points. One suggestion to find an approximate solutions to $\lambda = r_2(\lambda)$, but solving this takes some time, a bit too much for me. Hope this helps.