Solution (source code)

= Solution

With the <dimensionless variables> $x=r/r_s$ and $y=u/c_s$, <mass conservation> and part (c) give
$$
\frac\rho{\rho_0}=\frac{e^{3/2}}{x^2y}.
$$
Substitution into the <Bernoulli function>, together with $GM/r=2c_s^2/x$, eliminates every dimensional parameter and gives the <transcendental equation>
$$
\boxed{\frac{y^2}{2}-\ln y-2\ln x+\frac32-\frac2x=0}.
$$
Equivalently,
$$
\boxed{y^2-\ln(y^2)=\frac4x+4\ln x-3},
$$
or
$$
\boxed{y^2e^{-y^2}=x^{-4}e^{3-4/x}}.
$$
The <transonic branch> passes through $(x,y)=(1,1)$.