Plotting eigenvalues with different colours

Question

Bug introduced in version 11

---

I would like to plot all the eigenvalues (functions of Bz) of a matrix in the same plot, so as for Mathematica to understand that they are different functions and to plot in different colours.

Basically I used to do:

Plot[Evaluate[eigenvalues],{Bz, 0, 10}]

which used to work until I changed to Mathematica 11.0. Now it plots nothing.

Now I have to work around this by doing:

   Plot[eigenvalues /. Bz -> a, {a, 0, 10}]

which gives all the eigenvalues of the same colour. I can work with this but I'd like to understand why it's not plotting anymore.

Plotting the eigenvalues one by one I noticed that not all eigenvalues are not plottable, but only the ones where random '#' appeared in the expression.

For example, after running the code underneath, try to look at

eigenvalues[[7]]

It has #'s in it(why???), but if you do

 eigenvalues[[7]]/.Bz->100 //N

it still gives a finite answer!

--

FULL CODE:

MHz = 10^6;
h = 6.62607004*10^-34;  (* J\[Times]s *)
hbar = 1.054571800*10^-34; (* J\[Times]s *)

\[Mu]B = 9.27400968*10^-24; (* J\[Times]T^-1 *)
\[Mu]B = \[Mu]B/h*10^-4/10^6; (* G/MHz *)

splus[0] = {{0}} // SparseArray;
splus[S_] := 
 splus[S] = 
  SparseArray[
   Band[{1, 2}] -> 
    Table[Sqrt[S (S + 1) - M (M + 1)], {M, -S, S - 1, 1}], {2 S + 1, 
    2 S + 1}]
sminus[S_] := 
 Transpose[
  splus[S]] (* S- is just the transpose of S+ because have put the i \
in the definition of sx etc. *)
sx[S_] := (splus[S] + sminus[S])/2
sy[S_] := (splus[S] - sminus[S])/(2 I)
sz[S_] := SparseArray[
  Band[{1, 1}] -> Range[-S, S, 1], {2 S + 1, 2 S + 1}]
SparseIdentityMatrix[n_Integer /; n >= 1] := 
 SparseArray[Band[{1, 1}] -> 1, {n, n}]
id[S_] := id[S] = SparseIdentityMatrix[2 S + 1]

Ahfs = 6.093 ; (* MHz *)
Bhfs = 2.786 ; (* MHz *)
gJ = 4/3; (* ground *)
gI = -0.00014193489;
gS = 2.0023193043737;
gL = 0.99998627;

Ix = KroneckerProduct[sx[3/2], id[1/2], id[1]] ;(* tensor products *)


Iy = KroneckerProduct[sy[3/2], id[1/2], id[1]];
Iz = KroneckerProduct[sz[3/2], id[1/2], id[1]];
Itot = {Ix, Iy, Iz};

Sx = KroneckerProduct[id[3/2], sx[1/2], id[1]];
Sy = KroneckerProduct[id[3/2], sy[1/2], id[1]];
Sz = KroneckerProduct[id[3/2], sz[1/2], id[1]];
Stot = {Sx, Sy, Sz};

Lx = KroneckerProduct[id[3/2], id[1/2], sx[1]];
Ly = KroneckerProduct[id[3/2], id[1/2], sy[1]];
Lz = KroneckerProduct[id[3/2], id[1/2], sz[1]];
Ltot = {Lx, Ly, Lz};

Jtot = {Jx, Jy, Jz} = Stot + Ltot;
Ftot = {Fx, Fy, Fz} = Itot + Jtot;

Fs = Fx.Fx + Fy.Fy + Fz.Fz;
Js = Jx.Jx + Jy.Jy + Jz.Jz;

Hhfs = Ahfs*Sum[Itot[[i]].Jtot[[i]], {i, 1, 3}] + \[Mu]B*
    Bz*(gI*Iz + gS*Sz + gL*Lz); (* Bz in Gauss *)

{eigenvalues, eigenvectors} = 
  Eigensystem[
   SetPrecision[{Hhfs}, \[Infinity]]];(* Set that to \[Infinity] \
otherwise won't work *)

eigenvectors = Normalize /@ eigenvectors;

{eigenvalues, 
   eigenvectors} = {{eigenvalues[[1]], eigenvalues[[2]], 
    eigenvalues[[3]], eigenvalues[[5]], eigenvalues[[6]], 
    eigenvalues[[8]], eigenvalues[[9]], eigenvalues[[11]], 
    eigenvalues[[13]], eigenvalues[[14]], eigenvalues[[16]], 
    eigenvalues[[18]], eigenvalues[[19]], eigenvalues[[20]], 
    eigenvalues[[22]], eigenvalues[[24]]}, {eigenvectors[[1]], 
    eigenvectors[[2]], eigenvectors[[3]], eigenvectors[[5]], 
    eigenvectors[[6]], eigenvectors[[8]], eigenvectors[[9]], 
    eigenvectors[[11]], eigenvectors[[13]], eigenvectors[[14]], 
    eigenvectors[[16]], eigenvectors[[18]], eigenvectors[[19]], 
    eigenvectors[[20]], eigenvectors[[22]], eigenvectors[[24]]}};

Answer 1

At this time, this is only an extended comment.

When playing around with your problem, I noticed that the first 8 eigenvalues from the full eigenvalues list returned by Eigensystem plot just fine, with colors, using the expected syntax:

Plot[Evaluate@eigenvalues[[1 ;; 8]], {Bz, 0, 10}]

All the other ones don't plot at all:

MapIndexed[
  Plot[
    #1, {Bz, 0, 10},
    Epilog -> Inset[Style[First@#2, Red, 18], Scaled[{0.9, 0.4}]]
  ] &,
  eigenvalues
]

---

But there seems to be nothing inherently wrong with these expressions, nor does this behavior depend on the presence of Root objects: those are also present in the first 8, and those plot fine!

In fact, one can successfully reconstruct the full plot point by point:

ListPlot@
 Table[{Bz, eigenvalues[[i]]}, {i, Length@eigenvalues}, {Bz, 0, 10}]

---

After confirming that the original code works with no issues in MMA 10.4, I also found that Plot seem to calculate at least some of the points to display; it just doesn't display them:

Function[{arg},
 ListPlot @ #[[-1, -1]]& @
  Reap@
    Plot[
      arg, {Bz, 0, 10},
      EvaluationMonitor :> Sow[{Bz, arg}],
      PlotRange -> All
    ]
] /@ eigenvalues

Something funny is still going on, though: notice that some of the plots do not cover the whole requested range of Bz; a manual larger PlotPoints setting will force the calculation of all points, but nonetheless something seems off.