Solution (source code)

= Solution

Normalize all fluxes by $\pi$, rather than including that factor in only some of them. The specific plume <volume flux>, <momentum flux> and <buoyancy flux> are
$$
\boxed{Q=b^2w,\qquad M=b^2w^2,\qquad B=b^2wg'=B_s.}
$$
The ambient transport is $-Q$ and its specific <momentum flux> is
$$
M_a=(R^2-b^2)v^2=\frac{Q^2}{R^2-Q^2/M}.
$$
Put $s=b^2/R^2=Q^2/(R^2M)$, with $0<s<1$. The total specific <momentum flux> is
$$
J=M+M_a=\frac{M}{1-s},\qquad b=\frac{Q}{\sqrt M},\qquad w=\frac MQ,\qquad g'=\frac{B_s}{Q}.
$$
The first differential equation follows immediately from relative-speed <fluid entrainment>:
$$
\boxed{Q'=\frac{2\alpha\sqrt M}{1-s}.}
$$
For the second, differentiate $J=M/(1-s)$ and use $s'=2sQ'/Q-sM'/M$. Since $J'=b^2g'=B_sQ/M$, this yields
$$
\boxed{(1-2s)M'=\frac{B_sQ}{M}(1-s)^2-\frac{2Q}{R^2}Q'.}
$$
This is the requested two-equation system, with $B=B_s$ and zero total <volume flux> as the additional conservation conditions. Near the pure source, $Q,M\to0$ and $s\to0$. Eliminating height in that limit gives
$$
\frac{dM}{dQ}\sim\frac{B_sQ}{2\alpha M^{3/2}},\qquad M^{5/2}\sim\frac{5B_s}{8\alpha}Q^2.
$$
The zero <integration> constant selects a <pure plume>; a nonzero imposed source <momentum> would select a different branch and generally a different turning height. The leading radius is $b\sim6\alpha z/5$, which recovers the unconfined pure-plume limit near the source.