How to speed up this code?

I have Mathematica code for numerically evaluating integrals but it takes hours to run. Is there any way to improve its speed (I am concerned about numerical accuracy, so MonteCarloRule is not an option).

Here is the code:

Xi[R_, z_, m_, rc_] := rc^m/(rc^m + (R^2 + z^2)^(m/2))

Dz[R_, z_] := (R^2 + z^2)^(1/3);

Df[R_?NumericQ, z1_?NumericQ, z2_?NumericQ, m_?NumericQ,
rc_?NumericQ] :=
NIntegrate[(z1 - z) (Xi[0, z, m, rc] + Xi[R, z, m, rc]), {z, 0,
z1}] + NIntegrate[(z2 -
z) (Xi[0, z, m, rc] + Xi[R, z, m, rc]), {z, 0, z2}] -
NIntegrate[(Abs[z2 - z1] - z) Xi[R, z, m, rc], {z, 0,
Abs[z2 - z1]}]

DfD[R_?NumericQ, z1_?NumericQ, z2_?NumericQ, m_?NumericQ,
rc_?NumericQ] :=
NIntegrate[Xi[0, z, m, rc], {z, 0, z1}] +
NIntegrate[Xi[R, z1 - z, m, rc], {z, 0, z2}]

SDfD[R_?NumericQ, z1_?NumericQ, z2_?NumericQ, m_?NumericQ,
rc_?NumericQ] :=
NIntegrate[Xi[0, z, m, rc] + Xi[R, z2 - z, m, rc], {z, 0,
z1}] + NIntegrate[Xi[0, z, m, rc] + Xi[R, z1 - z, m,
rc], {z, 0, z2}]

XiI[R_?NumericQ, S_?NumericQ, m_?NumericQ,
rc_?NumericQ, Sigma_?NumericQ] :=
NIntegrate[
1/Sqrt[Dz[R,
z1 - z2]] Exp[-(z1 +
z2)] Exp[Sigma^2 Df[R, z1, z2, m,
rc]] (1 + Sigma^2 Xi[R, z1 - z2, m,
rc] - Sigma^2*
SDfD[R, z1, z2, m, rc] + Sigma^4 DfD[R, z1, z2, m,
rc] DfD[R, z2, z1, m, rc]), {z1, 0, S}, {z2, 0, S}]

At the end, I need to evaluate the following

XiI[0, 0.1, 1/3, 0.001, 0.01] -
ParallelTable[XiI[10^R, 0.1, 1/3, 0.001, 0.01], {R, -3, 0,
0.1}].

Unfortunately, this took about 12 hours to complete. Is there any trick to speed this up without compromising numerical accuracy?

• Have you tried Parallelize[]? By the way: I strongly suggest you avoid using upper-case-initial names for your variables and functions because they can conflict with Mathematica's internal names. – David G. Stork Dec 21 '16 at 2:47
• In my experience, there usually is more to be gained by improving the algorithm than by optimizing the code. For instance, if DfD[R, z1, z2, m, rc] is a slowly varying function of z1 and z2` over {0, S}. it may be advantageous to compute an array of values of DfD vs {z1, z2} and from it create an InterpolationFunction prior to performing the XiI integration. – bbgodfrey Dec 21 '16 at 5:26
• But the problem is DfD is not only function of {z1, z2} but also of R. – konstant Dec 21 '16 at 6:15