rmax=10;nlow=-500;nup=-1;eguess=-0.1;
epsilon=N[10^(-2)];tol=N[10^(-5)];

e[n_]:=eguess + n*epsilon;


(*************************************************************************)
l=1;iota=0.01;
(*sol[n_]:=NDSolve[{-1u'[r]-(r/2)u''[r]-(1-0.5*l*(l+1)/(r+0.0001))u[r]==r e[n]u[r],u[rmax]==rmax*Exp[-rmax],u[rmax-epsilon]==(rmax-epsilon)*Exp[-(rmax-epsilon)]},u,{r,epsilon/5,rmax}];*)
(*sol[n_]:=NDSolve[{-1u'[r]-(r/2)u''[r]-(1-0.5*l*(l+1)/(r+0.0001))u[r]==r e[n]u[r],u[rmax]==rmax*Exp[-rmax],u[rmax-epsilon]==(rmax-epsilon)*Exp[-(rmax-epsilon)]},u,{r,epsilon/5,rmax}];*)
sol[n_]:=NDSolve[{-r*u'[r]-((r^2)/2)u''[r]-(r-iota*0.5*l*(l+1))u[r]==r^2 e[n]u[r],u[rmax]==rmax*Exp[-rmax],u[rmax-epsilon]==(rmax-epsilon)*Exp[-(rmax-epsilon)]},u,{r,epsilon/5,rmax}];


(*************************************************************************)


usoldum[n_,r_]:=Evaluate[u[r]/.sol[n]];
usol[n_,r_]:=usoldum[n,rr][[1]]/.{rr->r};
usolzero[n_]:=usol[n,rr]/.{rr->0};
tusol[n_]:=Table[{r,usol[n,r]},{r,0,rmax/4,rmax/(4*20)}];
tusolzero[ndum1_,ndum2_]:=Table[{n,usolzero[n]},{n,ndum1,ndum2,(ndum2-ndum1)/30}];
n1=nlow;n2=nup;nave=N[(1/2)(n1+n2)];
Print["n1=",n1," usolzero[n1]=",usolzero[n1]];
Print["n2=",n2," usolzero[n2]=",usolzero[n2]];
ListPlot[tusolzero[nlow,nup],PlotLabel->"u[0] vs n"];

(*root finding*)
While[Or[usolzero[nave]>tol,(n2-n1)>tol],

If[usolzero[n1]*usolzero[nave]>0,n1=nave];
If[usolzero[n2]*usolzero[nave]>0,n2=nave];
nave=N[(1/2)(n1+n2)];

Print["n1=",n1," usolzero[n1]=",usolzero[n1]];
Print["n2=",n2," usolzero[n2]=",usolzero[n2]];
Print["nave=",nave," usolzero[nave]=",usolzero[nave]];
Print[ListPlot[tusolzero[n1,n2],PlotRange->{{nlow,nup},{-100,100}}]];
Print[" "];
];(*end while*)

Print["exit"];
Print["nave=",nave];
Print["e[nave]=",e[nave]];
Print["usolzero[nave]=",usolzero[nave]];
ListPlot[tusolzero[nave-5,nave+5]]
Plot[usol[nave,r],{r,0,rmax},PlotRange->{{0,rmax},All},PlotLabel-> {"u[r] vs r; e=",eguess + nave*epsilon}];
Manipulate[Plot[usol[n,r],{r,0,rmax},PlotRange->{-20,20},PlotLabel->{" ur[n] vs. r; n=",n,"e[n]=",e[n]}],{n,n1,n2,IntegerPart[(n2-n1)/20]}]

rmax=10;nlow=-500;nup=-1;eguess=-0.1;
epsilon=N[10^(-2)];tol=N[10^(-5)];

e[n_]:=eguess + n*epsilon;
(*************************************************************************)
l=1;iota=0.01;
(*sol[n_]:=NDSolve[{-1u'[r]-(r/2)u''[r]-(1-0.5*l*(l+1)/(r+0.0001))u[r]==r e[n]u[r],u[rmax]==rmax*Exp[-rmax],u[rmax-epsilon]==(rmax-epsilon)*Exp[-(rmax-epsilon)]},u,{r,epsilon/5,rmax}];*)
(*sol[n_]:=NDSolve[{-1u'[r]-(r/2)u''[r]-(1-0.5*l*(l+1)/(r+0.0001))u[r]==r e[n]u[r],u[rmax]==rmax*Exp[-rmax],u[rmax-epsilon]==(rmax-epsilon)*Exp[-(rmax-epsilon)]},u,{r,epsilon/5,rmax}];*)
sol[n_]:=NDSolve[{-r*u'[r]-((r^2)/2)u''[r]-(r-iota*0.5*l*(l+1))u[r]==r^2 e[n]u[r],u[rmax]==rmax*Exp[-rmax],u[rmax-epsilon]==(rmax-epsilon)*Exp[-(rmax-epsilon)]},u,{r,epsilon/5,rmax}];
Print[Plot[e[n],{n,nlow,nup}]];

e[nlow]

e[nup]

n=1;rmax=10;l=1;epsilon=0.01;
en=-0.1;iota=1;
sol:=u[r]/.DSolve[{-r*u'[r]-((r^2)/2)u''[r]-(r-iota*0.5*l*(l+1))u[r]==r^2 en*u[r],u[rmax]==rmax*Exp[-rmax],u[rmax-epsilon]==(rmax-epsilon)*Exp[-(rmax-epsilon)]},u[r],r][[1]];
f[rdum_]:=sol/.{r->rdum};

Print[f[r]]
(*
Plot[f[r],{r,0.01,10}]
*)
(*
For[i=1,i≤5,i++,
en=i*0.1;
Plot[f[r],{r,0.01,10}]
]
*)

f[0.0001]