Plot of spectrum of bound states of 2D well

(* Short notation for Bessel functions and their derivatives *)
j[n_,x_]:=BesselJ[n,x]
k[n_,x_]:=BesselK[n,x]
dj[n_,x_]:=(1/2) (BesselJ[-1+n,x]-BesselJ[1+n,x])
dk[n_,x_]:=(1/2) (-BesselK[-1+n,x]-BesselK[1+n,x])

(* style for different axial modes *)
style[0]:=Black
style[1]:=Red
style[2]:=Blue
style[3]:=Green;
style[4]:=Yellow
style[5]:=Orange
style[6]:=Purple

(* Define plot function *)
plot[n_]:=ContourPlot[{(Sqrt[1+ϵ]((dj[n,ξ0 Sqrt[1+ϵ]])/(j[n,ξ0 Sqrt[1+ϵ]]))==Sqrt[-ϵ]((dk[n,ξ0 Sqrt[-ϵ]])/(k[n,ξ0 Sqrt[-ϵ]])))},{ξ0,0,10},{ϵ,-1,0},Exclusions ->{j[n,ξ0 Sqrt[1+ϵ]]==0},PlotPoints->100,MaxRecursion->4,ContourStyle->style[n]]

(* Collect the plots *)
Show[{plot[0],plot[1],plot[2],plot[3],plot[4],plot[5],plot[6]},FrameLabel->{"(Subscript[ξ, 0])","ϵ"}]

(* Comment: the bound state for n=0 and small ξ0 is not found due to numericls*)