Coulomb Wave Function

Coul[n_, l_, z_, 
        r_]  =    Sqrt[
            z (n + l)!/(n - l - 1)!]/(n (2  l + 
                    1)!) (2 z r/n)^(l + 1) Exp[-z r/
              n] Hypergeometric1F1[l + 1 - n, 2 l + 2, 
          2 z r/n]

P10  =  Coul[1, 0, Z, r]

Do[
      Print[StringForm["P[``,``](r) = ``", n, l, Coul[n, l, Z, r]]], {l, 
        0, 4}, {n, l + 1, 5}]

Clear[Z]

Plot[{Coul[4, 0, 1, r], Coul[4, 1, 1, r], Coul[4, 2, 1, r], 
          Coul[4, 3, 1, r]}, {r, 0, 60}, 
        PlotStyle -> {RGBColor[0, 0, 0], RGBColor[1, 0, 0], 
            RGBColor[0, 1, 0], RGBColor[0, 0, 1]}, 
        PlotPoints -> 200];

Rsquare  =  
      Table[Integrate[r^2 Coul[4, l, Z, r]^2, {r, 0, Infinity}, 
          Assumptions -> Re[Z] > 0], {l, 0, 3}]

Table[8 (80 - 3 l (l + 1) + 1)/Z^2, {l, 0, 3}]

Rsquare  =  
      Table[Integrate[ Coul[4, l, Z, r]^2/r, {r, 0, Infinity}, 
          Assumptions -> Re[Z] > 0], {l, 0, 3}]

Table[
      Integrate[Coul[3, 2, Z, r] Coul[n, 2, Z, r], {r, 0, Infinity}, 
        Assumptions -> Re[Z] > 0], {n, 3, 10}]