The default gamma value has been upped to 1e4 which is enough to strongly constrain the solution to have zero slope at the core region. Further, the initial guess has been changed from a series expansion of theta to a simple quadratic that is one at origin and zero at the polytrope radius. This is faster to evaluate and seems to work just as well.