Why this code is considerably slower in version 9?

Question

I have a Mathematica code for computing and plotting the density which corresponds to a given potential. The following code was written few years ago using version 5.2 of the Mathematica. Recently, I tried this code in the latest version 9 and I observed something very disturbing. Using version 5.2 it needs about 5.5 sec in order to produce the contours. On the other hand, in version 9 I have to wait about a minute (ten times longer than in 5.2!) using of course the same PC. Why this is happening? Do I have to make any modification to make it more compatible with version 9?

Here is the code:

Clear["Global`*"];
Off[General::spell];

Vn = -(Mn/Sqrt[r^2 + z^2 + cn^2]);
Vd = -(Md/Sqrt[b^2 + r^2 + (a + Sqrt[h^2 + z^2])^2]);
Vh = v0^2/2*Log[r^2 + β*z^2 + ch^2];
Vt = Vn + Vd + Vh;
Vr = D[Vt, r];
Vrr = D[Vt, {r, 2}];
Vzz = D[Vt, {z, 2}];
ρ = 2.325/(4*π*100)*(Vrr + Vr/r + Vzz);
ρr = ρ /. {z -> 0};

Mn = 400; cn = 0.25;
Md = 5000; b = 6; a = 3; h = 0.2;
v0 = 20; β = 1.3; ch = 8.5;
ρ0 = ρr /. {r -> 8.5};

Print["-----------------"]
Print["ρ0 = ", ρ0]
Print["-----------------"]

S0 = ContourPlot[ρ, {r, -20, 20}, {z, -20, 20}, Contours -> 10, 
PlotPoints -> 200, ImageSize -> 550]

Answer 1

I think the PlotPoints needed in versions >=6 dropped for comparable quality, and so you might do better setting it around 60 or so.

Alternatively, keep it at 200 but set PerformanceGoal->"Speed"; I think that will be comparable to 5.2 behavior.