We're looking for solutions to
h m' m = 2 r m' - r^{2} m''
where h=(Gμ)/(RT) (G the universal gravitational constant, R the ideal gas constant, T the temperature of the (isothermal) gas, and μ the molecular weight); e.g. for air (μ=29g/mol) at 300K, h~8E-16m/kg. Since we want m to be the mass within radius r, physical solutions should have m>0 and m'>0.
First notice that there is one linear solution,
m = (2/h) r .
This gives a pressure
P = m'/A = 1/(2 pi) r^{-2} ,
diverging as r→0; so this is nonphysical for an all-gas sphere. (It may be reasonable if the boundary conditions have gas only for some r>r_{0}, e.g. the atmosphere of a planet, though real planet atmospheres are usually not isothermal.)
But let's write
m(r) = (2/h) r (1 + u(r)) ;
this gives us a differential equation for u,
0 = r^{2} u'' + 2 r (1 + u) u' + 2 (1 + u) u .
For physical solutions m>0 and so 1+u>0.
We can make this an autonomous differential equation by changing variables from r to t=log r: u'=(dt/dr)u'_{t}=(1/r)u'_{t} and u''=(dt/dr)(du'/dt)=1/r^{2}(u''_{tt}-u'_{t}) (decorating t derivatives with a subscript t to distinguish them from r derivatives), so
0 = u''_{tt} + (1 + 2u) u'_{t} + 2 (1 + u) u .
Since this is autonomous, it can be immediately converted to a first-order differential equation, though it's not obvious to me that this actually helps much.
More interesting to me is that this looks a lot like a damped harmonic-oscillator equation, but with varying coefficients: undamped (squared) frequency 2(1+u) and damping (1+2u). Since both of these are positive, this equation should be stable with increasing t, i.e. it should only support solutions with u→0 as t→∞. That is, if we pretend the coefficients 2(1+u) and (1+2u) are nearly constant, the solutions to the resulting equation are linear combinations
u ~ A exp(σ_{1}t) + B exp(σ_{2}t)
with coefficients
σ_{1,2} = -(1/2) ( (1+2u) ± sqrt(-7 - 4u + 4u^{2}) ) .
For small |u| these are complex conjugates, giving an underdamped solution decaying roughly as exp(-t/2)*sin(t*sqrt(7)/2+φ). This means that for sufficiently large t (hence r), solutions should converge toward u→0 and hence toward m ~ (2/h) r.
The solutions I think you're most interested in are those which are gas all the way to r=0. Near r=0 for a gas sphere the pressure should be nearly constant (since the gravitational gradient is small), so m=ρr^{3} for some near-center density rho.
So then we expect the full solution to look roughly like
m ~ min (ρr^{3} , (2/h) r)
with some oscillations in the transition region.