Solving a second order ODE with ParametricNDSolveValue

I am trying to numerically solve a nonlinear ODE on the interval 0 to R, using ParametricNDSolveValue. The boundary conditions are : x' == 0 at r == R, x' == 0 at r == 0, and x at R==0 .c0 is a parameter.

my function:

x''[r] + x'[r]/r == -c0 Exp[-x[r]]+357/25 for 0 <= r <= a,
x''[r] + x'[r]/r == -c0 Exp[-x[r]]        for a < r <= R

Here is my code that I just did it:

R = 78958914/46978895;

q = 100;
f[r_] := Piecewise[{{357/2500*q, 0 <= r <= a}, {0, a < r <= R}}]
sol = ParametricNDSolveValue[{x''[r] + x'[r]/r == -(357 \[Pi])/1250 c0 Exp[-x[r]] +
f[r], x[\[Epsilon]] == x0, x'[\[Epsilon]] == 0,
WhenEvent[r == 1, x'[r] -> x'[r]]}, {x, x'}, {r, \[Epsilon],
R}, {c0}, Method -> "StiffnessSwitching",
WorkingPrecision -> 30]; ff =
FindRoot[Last[sol[x0]][R], {x0, -11, 0}, Evaluated -> False][[1, 2]]

It gives my a solution, but it is not correct, why?, the area under the curve should equal to q=100 and when I integrated over it, it gave me 3.93993*10^-8! by using q = 2 \[Pi] ff NIntegrate[r Exp[-First[sol[ff]][r]], {r, 0, R}]. • Comments are not for extended discussion; this conversation has been moved to chat. – Kuba Jun 11 at 9:12