# How to solve this ODE by NDSolve?

I am trying to solve an ODE using NDSolve, but I get a division by zero error. I am not sure how to fix this and obtain a solution:

n=3;
g[e_,r_]=e*(2r^3/3 -r^2);
c=1/2;
H[e_,r_]=Integrate[s^2*(c-2*g[e,s]),{s,0,r}];
odeN= {(n-2)*r*f'[r]*Cot[f[r]]/Sqrt[1+(r*f'[r])^2]-D[r^3*f'[r]/Sqrt[1+(r*f'[r])^2],r]-(n-3)f'[r]/Sqrt[1+(r*f'[r])^2]==H[1,r], f[0] == 1,f'[0]==0}
solN= NDSolve[odeN, f, {r, 0, 20}];
pN=Plot[f[r]/. solN, {r, 0, 20}, PlotStyle -> {Blue}];
Show[pN,ImageSize->Large]


Is there a way around this issue?

• Just for the test, change the initial point to a sufficiently small number instead of $0$, e.g. f[0.0001] == 1, f'[0.0001] == 0 and solN = NDSolve[odeN, f, {r, 0.0001, 20}] then check the solution. It appears to have a singularity at 0 Commented Jan 24, 2023 at 18:57
• It works but I get a warning: singularity or stiff system expected. Are there methods that deal with stiff/singular systems? Commented Jan 24, 2023 at 19:15

ClearAll["Global*"]
ic = {f[$$MachineEpsilon] == 1, f'[$$MachineEpsilon] == 0}
`