NDSolve divide-by-zero trouble

Determining n'[r] from

Solve[eulerEq, Derivative[1][n][r]][[1, 1, 2]]
(* -((2 n[r] (-1 + r^3 n[r]^2))/(r (-1 + r^4 n[r]^(12/5)))) *)

we see that it is indeterminate at r = 1, and it is for this reason that NDSolve fails there. Nonetheless, the limiting value can be obtained as follows,

Limit[% /. n[r] -> 1 + a (r - 1), r -> 1]
(* -((5 (3 + 2 a))/(10 + 6 a)) *)
sl = Solve[% == a]
(* {{a -> 1/6 (-10 - Sqrt[10])}, {a -> 1/6 (-10 + Sqrt[10])}} *)

and used to define boundary conditions very near r = 1.

dr = .0001;
s1 = NDSolveValue[{eulerEq, n[1 + dr] == 1 + dr a /. First[sl]}, n, {r, 1 + dr, 10}];
s2 = NDSolveValue[{eulerEq, n[1 + dr] == 1 + dr a /. Last[sl]}, n, {r, 1 + dr, 10}];
s3 = NDSolveValue[{eulerEq, n[1 - dr] == 1 - dr a /. First[sl]}, n, {r, .5, 1 - dr}];
s4 = NDSolveValue[{eulerEq, n[1 - dr] == 1 - dr a /. Last[sl]}, n, {r, .5, 1 - dr}];

p1 = Plot[{s1[r], s2[r]}, {r, 1 + dr, 10}, PlotRange -> {{0, 10}, {0, 5}}];
p3 = Plot[{s3[r], s4[r]}, {r, .5, 1 - dr}, PlotRange -> {{0, 10}, {0, 5}}];
Show[p1, p3, AxesLabel -> {r, n}]

enter image description here


bbgodfrey has already given a good answer, here's just a not that rigorous but cheap and accurate enough way to resolve the problem:

S = NDSolveValue[{eulerEq, n@1 == 1 + #}, n[r], {r, 0.5, 10}, 
    Method -> "StiffnessSwitching"] & /@ ({1, -1} $MachineEpsilon)

Plot[S, {r, 0.5, 10}, PlotRange -> All]

enter image description here


Using the limiting values from bbgodfrey's answer, we can take a more direct approach with NDSolve using "EquationSimplification" -> "Residual":

sol = Join @@@
       (NDSolve[{eulerEq, initcond, n'[1] == #}, 
         n[r], {r, 0.5, 10},
         Method -> {"EquationSimplification" -> "Residual"}] & /@
          {1/6 (-10 + Sqrt[10]), 1/6 (-10 - Sqrt[10])});

Plot[n[r] /. sol // Evaluate, {r, 0.5, 10}, PlotRange -> All, 
 AxesOrigin -> {0, 0}]

Mathematica graphics