L  =  1
J  =  3/2
m  =  206.771
Z  =  82
Rad  =  0.0001343
Energy  =  4.5362

ϵ  =  Energy*(1000000/27.1139565)

(*Defining the exterior reduced radial wavefunction *) 
W[
        z_]  =  (E^(- z / 2)) (z^(L  +  
              1)) HypergeometricU[1 + L  -  η, 2  L  +  2, 
          z]

    η  =  Z (Sqrt[(m/(2 ϵ))])

    z  =  (Sqrt[8 m ϵ]) r

Rexterior[r_]  =  W[z]

(*Defining the interior reduced radial wavefunction *) 
Rinterior[
        r_]  =  (r^(L  +  
              1)) (E^(- 
                4 ϵ m r^2 β^2)) 
        Hypergeometric1F1[-γ,  L  +  3/2,  
          8 ϵ m β^2 r^2]

    β  =  ((((Z/(64 m Rad^3 ϵ^2))))^0.25)

    Λ  =  ((3 Z )/(8 Rad ϵ)) -  (1/4)

    γ  =  (1/2) ( (Λ/(2 β^2)) -  
            (1/2) -  (L  +  1))

Const  =  
      Rinterior[Rad]/Rexterior[Rad]

R[r_]  =  If[r < Rad,  Rinterior[r], Rexterior[r] Const]

Norm  =  
      NIntegrate[R[r]^2,  {r, 0, ∞}]

    R[r_]  =  R[r]/(Sqrt[Norm])

V[r_]  =  
      If[r < Rad,  - (Z/Rad) ((3/2)  -  
              r^(2/(2 Rad^2))),  - (Z/r)]

Vbar  =  
      NIntegrate[R[r]^2 V[r],  {r, 0, ∞}]

    Vsquare  =  
      NIntegrate[(R[r] V[r])^2,  {r, 
          0, ∞}]

    Darwin  =  
      NIntegrate[
        R[r]^2 ((3 Z )/Rad)^3,  { r, 0, 
          Rad}]

    ξ[r_]  =  
      If[L == 0, 0, If[r > Rad, (Z/r)^3,  (Z/Rad)^3]]

    SO  =  NIntegrate[R[r]^2 ξ[r], {r, 0, ∞}]

(*Shifting the energy*) 
α  =  0.007297353013

    Efs  =  - ϵ  -  α^(2/(2 m)) (ϵ^2  +  2 ϵ Vbar  +  
              Vsquare)  +  α^(2/(8 m^2)) 
          Darwin  +  α^(2/(4 m^2)) 
          SO (J (J + 1)  -  L (L  +  1)  -  
              (3/4))

    EfsMev  =  Efs*27.21139565/1000000