Using NMinimize with Table to plot an energy

I would like to reproduce the solid blue line of Figure 2a in this article. For this, I have used the code found here. Notice, however, that I am not succeeding. Could you help me? Here is my code:

l[wx_, wz_] := wz/wx;

f[X_] := (1 + 2 X^2)/(1 - X^2) -
3 X^2 ArcTanh[Sqrt[1 - X^2]]/(1 - X^2)^(3/2);

32*4 Pi*h^2*as/(3 M) Sqrt[as^3/Pi] (1 + 3 add^2/(2 as^2));

az[h_, M_, wz_] := Sqrt[h/(M*wz)];

etot[x_, y_, as_] :=
2 az[h, M, wz]^2/x^2 + az[h, M, wz]^2/y^2 +
x^2/(8 l[wx, wz]^2 az[h, M, wz]^2) + y^2/(16 az[h, M, wz]^2) +
8 n*az[h, M, wz]^3/(Sqrt[2 Pi] x^2 y) (as/az[h, M, wz] -
128 n^(3/2) yqf[h, M, as,
add]/(25 Sqrt[5] Pi^(9/4) x^3 y^(3/2) h*wz)

a0 = 5.29*10^(-11);
wx = 2 Pi*45;
wz = 2 Pi*133;
h = 1.054*10^(-34);
M = 2.72*10^(-23);
n = 15*10^3;

Clear[x, y, as]
ddata = Table[{x, y, as} /.
Last[NMinimize[{etot[x, y, as], x > 0 && y > 0}, {x, y}]], {as,
40*a0, 90*a0, 10*a0}]

mminima = Re[etot[#[[1]], #[[2]], #[[3]]]] & /@ ddata

uvalues = Table[as, {as, 40*a0, 90*a0, 10*a0}];
ListPlot[Transpose[{uvalues/a0, mminima}], Joined -> True,
Frame -> True]


Here, $$M$$ is Dy-164 mass and $$h$$ is the $$\hbar$$.

• It seems to be working. What is the problem that you observed? – MelaGo Oct 8 '19 at 0:49
• @MelaGo My plot has not negative values for instance. – Dinesh Shankar Oct 8 '19 at 2:15
• Are you sure M shouldn't be 2.72*10^-25? – MelaGo Oct 8 '19 at 3:50

I think the error was the value for M. (Also cleaned up the code a bit.)

ClearAll["Global*"]

a0 = 5.29*10^(-11);
wx = 2 Pi*45;
wz = 2 Pi*133;
h = 1.054*10^(-34);
(*  M=2.72*10^(-23);  *)
M = 163.9 1.66 10^-27;
n = 15*10^3;
az = Sqrt[h/(M*wz)];
l = wz/wx;

f[X_] := (1 + 2 X^2)/(1 - X^2) - 3 X^2 ArcTanh[Sqrt[1 - X^2]]/(1 - X^2)^(3/2);

yqf[as_] := 32*4 Pi*h^2*as/(3 M) Sqrt[as^3/Pi] (1 + 3 add^2/(2 as^2));

etot[x_, y_, as_] := 2 az^2/x^2 + az^2/y^2 + x^2/(8 l^2 az^2) + y^2/(16 az^2) +
8 n*az^3/(Sqrt[2 Pi] x^2 y) (as/az - add*f[x/y]/az) +
128 n^(3/2) yqf[as]/(25 Sqrt[5] Pi^(9/4) x^3 y^(3/2) h*wz)

ddata = Table[
minsol = NMinimize[{etot[x, y, as ], x > 0 && y > 0}, {x, y}];
Re[{x, y, as/a0 , minsol[[1]]} /. Last[minsol]], {as, 40 a0, 90 a0, a0}]

(* {{2.54948*10^-7, 5.21789*10^-6, 40., -87.6828},
{2.61239*10^-7, 5.25204*10^-6, 41., -81.3895},
{2.67699*10^-7, 5.28499*10^-6, 42., -75.5096},
{2.74338*10^-7, 5.31675*10^-6, 43., -70.0129}, etc *)

ListPlot[ddata[[All, {3, 4}]], Joined -> True, PlotRange -> {Automatic, 10}]
`

This curve looks funny, but it does resemble the Max of the blue and black solid curves in Fig 2 in the paper.