Need help with "Filling" in a plot, generated with a table. (Wavefunctions)

Question

I'm trying to get a result like this:

f[n_, x_] := Abs[((1/Pi)^(1/4) HermiteH[n, x])/(E^(x^2/2) Sqrt[2^n n!])]^2

Plot[Evaluate@Append[Table[f[n, x] + n + 1/2, {n, 0, 7}], x^2/2], {x, -4, 4}, Filling -> gt; Table[n -> gt; n - 1/2, {n, 1, 8}]]

Which gives this.

Only using this code:

a[n_, mu_, delta_] := 1/2*(mu/delta^2*1/(1 + n) - n - 1)

psi[n_, r_, mu_, delta_] := (-1)^n*Pochhammer[2, n]*Exp[-delta*a[n, mu, delta]*r]*(1 - Exp[-delta*r])*Hypergeometric2F1[-n, n + 2*a[n, mu, delta] + 2, 2, 1 - Exp[-delta*r]]

En[n_, mu_, delta_] := -(mu/(2 delta (n + 1)) - 1/2 delta (n + 1))^2

V[r_, mu_, delta_] := -mu*Exp[-delta*r]/(1 - Exp[-delta*r])

cns[n_, mu_, delta_] := 1/Sqrt[(Pochhammer[2, n])^2/(delta*(a[n, mu, delta] + 1)*(2*a[n, mu, delta] + 1)*(2*a[n, mu, delta]))*Sum[Pochhammer[-n, k]*Pochhammer[n + 2*a[n, mu, delta] + 2, k]*Pochhammer[3, k]/(Pochhammer[2, k]*Pochhammer[2*a[n, mu, delta] + 3, k]*k!)*HypergeometricPFQ[{-n, n + 2*a[n, mu, delta] + 2, k + 3}, {2, 2*a[n, mu, delta] + k + 3}, 1], {k, 0, n}]]

npsi[n_, r_, mu_, delta_] := (-1)^n*Pochhammer[2, n]*cns[n, mu, delta]*Exp[-delta*a[n, mu, delta]*r]*(1 - Exp[-delta*r])*Hypergeometric2F1[-n, n + 2*a[n, mu, delta] + 2, 2, 1 - Exp[-delta*r]]



With[{mu = 6^2, delta = 1}, {Plot[Evaluate@Append[Table[npsi[n, r, mu, delta] - (mu/(2 delta (n + 1)) - 1/(2 delta (n + 1)))^2, {n, 1, 4}], V[r, mu, delta]], {r, 0, 5}, PlotRange -> gt; {-80, 1}, Filling -> gt; Axis]}]

(I hope I included all relevant code, let me know if I didn't).

Which gives this at the moment.

I realize I have the filling set to "Axis" but I can't seem to figure out what to put in for it and when I try it keeps giving me an error saying it's not a valid filling interval.

I'm looking for it to fill basically at the energy level for the given wavefunction to the peaks of the wavefunction it's plotting after.

Answer 1

something like this?

functionsToPlot[r_, mu_, delta_] := Table[
  N[npsi[n, r, mu, 
     delta] - (mu/(2 delta (n + 1)) - 1/(2 delta (n + 1)))^2],
  {n, 1, 4}
  ]

With[{mu = 6^2, delta = 1},
 Plot[
  {
   V[r, mu, delta],
   Evaluate[functionsToPlot[r, mu, delta]]
   },
  {r, 0, 5},
  PlotRange -> {-82, 1},
  Filling -> {
    2 -> -79,
    3 -> -36,
    4 -> -22,
    5 -> -14
    },
  ImageSize -> Large
  ]
 ]

I defined functionsToPlot just to make the code more clear to read. As a side note, you should avoid variable names starting with an uppercase letter, to prevent unexpected conflicts with built-in symbols.