One more tidbits to it. When I tried to plot it without Evaluate, it grabbed 20GB as "Cached Files". So, it was doing some heavy duty work on my SSD - although still there was plenty of real memory available -, but I realized it only after the facts. I do not know if it is still one of the side effects of the Frontend is being 32-bit. There were times a few years ago when I pushed the limits of this machine into 1% idle and constant swapping in and out of memory, virtual memory were going into hundreds of GBs, but even at those times the "Cached Files" never grew up to 20GB, see picture attached.

I would have thought that with immediate definitions (=) instead of delayed definition (:=) for the quantities rbmmsp, rbmmf would be equivalent to Evaluate, because the very complicated formulas are pre-calculated, before the insertion of numeric values. Surprisingly to me, there is no speed-up.

The lesson here seems to be like this: if we are plotting a function with a very lengthy algebraic formula, we had better feed it to Plot etc in expanded form, either with Evaluate or With (or, conceptually, copy-paste of the formula itself).

Thanks to those who brought this to our attention. In investigating the slowdown in recent versions of Mathematica, it looks like the slow evaluation is a result of VectorPlot3D needing to evaluate the nqvrbmmf[mu0, a, b, c, x, y, z] function at all the points that the plot samples. You can see that individual calls of nqvrbmmf are relatively slow:

In[5]:= Do[
  nqvrbmmf[1000000, 0.01`, 0.01`, 0.01`, x, y, 
    z] /. {x -> RandomReal[{-0.015`, 0.059`}], 
    y -> RandomReal[{-0.02`, 0.02`}], 
    z -> RandomReal[{-0.035`, 0.115`}]}, 200] // AbsoluteTiming

Out[5]= {19.7192, Null}

This slow evaluation time can be addressed by forcing evaluation of VectorPlot3D's first argument with Evaluate prior to the actual evaluation of the plot:

In[6]:= AbsoluteTiming[
 VectorPlot3D[
     Evaluate[nqvrbmmf[mu0, a, b, c, x, y, z]] , {x, -1.5*a, 
      5.5*a + 4*gap}, {y, -2*b, 2*b}, {z, -3.5*c, 1.5*c + R}, 
     VectorPoints -> {40, 20, 30}, BoxRatios -> Automatic, 
     Axes -> True, PlotLegends -> Automatic, 
     PlotRange -> {{-0.015, 0.059}, {-0.02, 0.02}, {-0.035, 0.115}}] //
     ReplaceAll[#, {mu0 -> 10^6, a -> 0.01, b -> 0.01, c -> 0.01, 
       gap -> 0.001, R -> 0.1}] & // Quiet;]

Out[6]= {1.43917, Null}

The problem here was actually in Mathematica 11.0.1, where it looks like VectorPlot3D's first argument was evaluating before it should, leading to the speedup and causing the confusion.

Just to be sure I submitted too. Thanks a lot. János

Thanks Marco. Then it is not just my nightmare. I am wondering if I should sent it to Support too.

Have you tried comparing the formulas for nqvrbmmf in the two versions? Are they the same formula?