ClearAll["Global`*"]

hbar=1.054*10^-34;
m0=1.67*10^-27;


ρ1=16000*10^-18;                                         
ρ2=2200*10^-18;
mass=1.16*10^-10;
R0=(mass/ρ1)^(1/3); 


Do[
R=R0 * 2^rr;
                                                                                                                                 
rc=0.1;
IR=8 (π^3) (R^2) ((1/2)-(1/2) (E^(-((R^2)/(2 (rc^2))))) (BesselI[0,((R^2)/(2 (rc^2)))]+BesselI[1,((R^2)/(2 (rc^2)))]));
Sffexp=2*10^-36;     

Do[
ϵ=2^nϵ;
Nlay=2^n;
c1=mass/((((Nlay+1)*ρ1+Nlay*ρ2))*π*R^2);                                                                                     
a=mass/(π R^2((Nlay+1)ρ1+Nlay ϵ ρ2));
b=ϵ a;
L=(Nlay+1)a+Nlay b;                        

ρkz=2/(kz Sin[kz(a+b)/2])(ρ1 Sin[kz a/2]Sin[kz (L+b)/2]+ρ2 Sin[kz b/2]Sin[kz (L-a)/2]);

Kz=kz^2*ρkz^2*Exp[-kz^2*rc^2];
IR=8 (π^3) (R^2) ((1/2)-(1/2) (E^(-((R^2)/(2 (rc^2))))) (BesselI[0,((R^2)/(2 (rc^2)))]+BesselI[1,((R^2)/(2 (rc^2)))]));
Iz=NIntegrate[Kz,{kz,-Infinity,+Infinity},MaxRecursion->40,Method->{GlobalAdaptive,MaxErrorIncreases->10000}];
dati=rc^3*IR*Iz/(π^(3/2)m0^2);

ULdati=Sffexp/(dati*hbar^2*2)*10^-12;

count+=1;
If[ULdati<best,best=ULdati;countbest=count;Print["Best at =",count];];
AppendTo[XY,{rr,2^n,nϵ,ULdati}];
,{nϵ,-2,2,2},{n,0,12,2}];
,{rr,3,5,1}]