a0=1000000; 
a1=-10000000;
a2=-5000000;
(*d0=0.005;*)
h=0.001;
rmin= 1;
rmax= 30000;



bd1[x_,t_]:=1+(-9/10)t x^(-2/3)+(81/280)t^2 x^(-4/3)+((1/6)-(27t^3/560))x^(-2)+((-27t/110)+(9 27t^4/(88 560)))x^(-8/3);
bd2[x_,t_]:=(-2/3)(-9/10)t x^(-5/3)+(-4/3)(81/280)t^2 x^(-7/3)+(-2)((1/6)-(27t^3/560))x^(-3)+(-8/3)((-27t/110)+(9 27t^4/(88 560)))x^(-11/3);


(*oklist={};*)
Do[y1[j]=NDSolve[{y''[x]+(-(1+x^2)^(-2)+h j(1+x^2)^(-4/3))y[x]==0,y[a0]==bd1[a0,h j],y'[a0]==bd2[a0,h j]},y,{x,a1,a0}]; 
y2[j]=y[a1]/.y1[j];
y3[j]=y[a2]/.y1[j];
y4[j]=y3[j][[1]]-y2[j][[1]];
(*flag[j]=If[-d0<y4[j]<d0,AppendTo[oklist,h j];"OK","NG"]; *)
(*Print["lambda=",h j,", ",flag[j],", the difference is ",y4[j]]*),{j,rmin,rmax}]
(*Print["candidates of eigenvalues are ",oklist]*)
ListPlot[Table[{h j,y4[j]},{j,rmin,rmax}]]