Psi[ξ_,η_,ϕ_,R_,α_]=Exp[-α*ξ*R/2]Cosh[α*η*R/2];

Vol[ξ_,η_,ϕ_,R_]=R^3/8(ξ^2-η^2);
V[ξ_,η_,R_]=-2/(R(ξ+η))-2/(R(ξ-η))+1/R;

TProlate[ξ_,η_,ϕ_,R_,α_]=-2/(R^2(ξ^2-η^2))*(D[(ξ^2-1)*D[Psi[ξ,η,ϕ,R,α],{ξ,1}],{ξ,1}]+D[(1-η^2)*D[Psi[ξ,η,ϕ,R,α],{η,1}],{η,1}]);

Tii[R_,α_]=Assuming[R>0&&α>0,Integrate[Vol[ξ,η,ϕ,R]*Psi[ξ,η,ϕ,R,α]*TProlate[ξ,η,ϕ,R,α],{ξ,1,∞},{η,-1,1},{ϕ,0,2π}]]
Vii[R_,α_]=Assuming[R>0&&α>0,Integrate[Vol[ξ,η,ϕ,R]*Psi[ξ,η,ϕ,R,α]*V[ξ,η,R]*Psi[ξ,η,ϕ,R,α],{ξ,1,∞},{η,-1,1},{ϕ,0,2π}]]
Sii[R_,α_]=Assuming[R>0&&α>0,Integrate[Vol[ξ,η,ϕ,R]*Psi[ξ,η,ϕ,R,α]*Psi[ξ,η,ϕ,R,α],{ξ,1,∞},{η,-1,1},{ϕ,0,2π}]]
En[R_,α_]=FullSimplify[(Tii[R,α]+Vii[R,α])/Sii[R,α]]
DEn[R_,α_]=FullSimplify[D[En[R,α],{α,1}]]

WykOptR=Array[0.&,{500,2}];
EnOptR=Array[0.&,{500,2}];
For[i=0,i≤499,i++,
R=0.001+i*0.1;
WykOptR[[i+1,1]]=R;
EnOptR[[i+1,1]]=R;
WykOptR[[i+1,2]]=FindRoot[DEn[R,α]==0,{α,1}][[1,2]];
EnOptR[[i+1,2]]=En[R,WykOptR[[i+1,2]]];
]
Clear[R];

WykOptFunkcR=Interpolation[WykOptR,Method->"Spline"];
EnOptFunkcR=Interpolation[EnOptR,Method->"Spline"];

Limit[En[R,α]-1/R,R->0]

Plot[{EnOptFunkcR[R],En[R,1]},{R,0.6,20},PlotRange->Full]

Plot[WykOptFunkcR[R],{R,3.5,25},PlotRange->Full]

αMin=FindMinimum[En[r,1],{r,1.2}]
αOptMin=FindMinimum[EnOptFunkcR[r],{r,1.2}]

(*Bonus: Rysujemy funkcje falowe. Wspolczynniki przejscia Prolate->Cartesian. Polozylem wszedzie R=1*)
ξc[x_,y_,z_]=Re[ArcCosh[I (Sqrt[(x^2)+(y^2)])+z]];
ηc[x_,y_,z_]=Im[ArcCosh[I (Sqrt[(x^2)+(y^2)])+z]];
ϕc[x_,y_,z_]=ArcTan[x,y];

PsiCartesian[x_,y_,z_,α_]=Exp[-α*ξc[x,y,z]/2]Cosh[α*ηc[x,y,z]/2];

Manipulate[Plot[PsiCartesian[0,0,z,α]^2,{z,-6,6},PlotRange->Full],{α,0.01,2.5}]

Manipulate[Plot3D[PsiCartesian[0,y,z,α]^2,{y,-10,10},{z,-10,10},PlotRange->Full],{α,0.1,3}]