Thomas-Fermi for \(TraditionalForm\`\(Al\^+\)\)

Following the discusions on the positive ions in the lecture notes, we need to find χ'(0) so that (Subscript[x, 0]) χ' ((Subscript[x, 0])) = -(Z - N)/
            Z = -1/13 for \(TraditionalForm\`\(Al\^+\)\).  After some trial-and-error, I find χ' (0) = -1.5881063075 works:

chisol = 
      NDSolve[{(1 + (4/3) (x^(3/2))) y''[x] + 
              4  (x^(1/2)) y'[x] + 
              If[x == 0, 
                0, y[x]/(x^(1/2)) (1 - (((1 + (4/3) 
                                  (x^(3/2))))^(3/2)) 
                        Max[y[x], 0]^(1/2))] == 0, 
          y[0] == 1, y'[0] == -1.5881062877}, y, {x, 0, 300}]; χ[
        x_] := (1 + 4/3*(x^(3/2))) y[x]; 
    Plot[χ[x] /. chisol, {x, 0, 300}]

Plot[χ[x] /. chisol, {x, 0, 30}]

Look for where χ((Subscript[x, 0])) = 0,

(Subscript[x, 0]) = 
      x /. FindRoot[Evaluate[χ[x] /. chisol], {x, 0, 100}]

Now make sure (Subscript[x, 0]) χ' ((Subscript[x, 0])) = -1/13,

Evaluate[Z x χ'[x] /. 
            chisol /. {x -> (Subscript[x, 0])} /. {Z -> 13}]

Using the definition e Φ(r) = ((Z e)/r) χ(
            r) and variables r =  (Z^(-1/3)) b x, b = (1/2) (((((3  π)/4)))^(2/
                3)) ((ℏ^2)/(m (e^2))), we define (in the atomic unit e = m = ℏ = 1)

Φ[
        r_] := -(Z/r) χ[((((Z^(-1/
                                3)) (1/2) (((((3  π)/4)))^(2/3))))^(-1)) 
                  r] /. {(Subscript[r, 0]) -> (Z^(-1/
                        3)) (1/2) (((((3  π)/4)))^(2/3)) 
                  (Subscript[x, 0])} /. chisol[[1]] /. {Z -> 13}

Plot[r Φ[r], {r, 0, 20}]

The Fermi energy is just (Subscript[ϵ, F]) = ((-e)/(Subscript[r, 0])), where (Subscript[r, 0]) = (Z^(-1/3)) 
            b (Subscript[x, 0]) = (Z^(-1/
                  3)) (1/2) (((((3  π)/4)))^(2/3)) 
            (Subscript[x, 0]) in the atomic unit.

ϕ[
        r_] := -(Z/r) χ[((((Z^(-1/
                                  3)) (1/2) (((((3  π)/4)))^(2/3))))^(-1)) r] - 
              (1/(Subscript[r, 0])) /. {(Subscript[r, 0]) -> (Z^(-1/
                        3)) (1/2) (((((3  π)/4)))^(2/3)) 
                  (Subscript[x, 0])} /. chisol[[1]] /. {Z -> 13}

(Z^(-1/3)) (1/2) (((((3  π)/4)))^(2/3)) 
        (Subscript[x, 0]) /. {Z -> 13}

Plot[r ϕ[r], {r, 0, 10}]

Indeed, ϕ(r) becomes (-1) ((e^2)/r) due to the overall positive charge beyond r > (Subscript[r, 0]) where the electron number density vanishes.  The coefficient is consistent with -1 within the tiny numerical errors.  To make this point clear, let us plot it on bigger scales,

Plot[r ϕ[r], {r, 0, 10}, PlotRange -> {-14, 0}, 
      PlotStyle -> RGBColor[0, 1, 0]]

Plot[r ϕ[r], {r, 0, 10}, 
      PlotRange -> {-1.01, -0.99}]