(*Finds the decreasing positive solutions of radial Neumann problem 
(Subscript[u, rr])+((n-1)/r) (Subscript[u, r]) + λ f(u)=0 for r<1, u'(0)=u'(1)=0, in n dimensions. Draws the bifurcation curve (λ, u(0))*)
Clear["`*"]
f[u_]=u (u-1) (7-u);
n=3;
u00=0;
delu=0.1;
nsteps=70;
tend= 1000;
h=$MachineEpsilon;
neumann={};
For[k=1,k ≤ nsteps,k++,
 u0=u00+k*delu;
sol=NDSolve[{x''[t]==-(n-1)*x'[t]/t-f[x[t]],x[h]==u0,x'[h]==0,WhenEvent[{x[t]*x'[t]==0},If[Abs[x'[t]]<10^(-8),AppendTo[neumann,{t^2,u0}]];"StopIntegration"]},x,{t,h,tend}];
]

ListPlot[neumann,Joined -> True,  AxesLabel -> {"λ","u(0)"}, PlotRange -> All, PlotStyle -> Thickness[0.008],PlotLabel->"Bifurcation Diagram for Neumann problem"]