Speeding up a slow indefinite integral

Question

I have found an alternative way to express an old solution:

PwD = (1/(4*Sqrt[Pi])) Integrate[(Exp[-(0.25/t)]/t^1.5)*
    Sum[Exp[-((n^2*Pi^2*t)/40000)] 
        Integrate[
            Exp[-((z^2*Tan[45*Degree]^2 + (200*n - z)^2)/(4*t))] + 
            Exp[-((z^2*Tan[45*Degree]^2 + (100 + 200*n - z)^2)/200)], 
            {z, -50, Plus[50]}
        ], {n, -4, 4}
    ], {t, 0, y}  
];

Solving a partial solution, I have:

PwD = (Exp[-(0.25/t)]/t^1.5)*Sum[Exp[-((n^2*Pi^2*t)/40000)]*
    Integrate[
        Exp[-((z^2*Tan[45*Degree]^2 + (200*n - z)^2)/(4*t))] + 
        Exp[-((z^2*Tan[45*Degree]^2 + (100 + 200*n - z)^2)/200)], 
        {z, -50, Plus[50]}
    ], {n, -4, 4}
];

Now, for plotting the indefinite integral in time (t) I had used the following codes:

function[y_] := (1/(4 Sqrt[π]))*NIntegrate[PwD, {t, 0.01, y}];
tabulate = Table[{y, function[y]}, {y, 0, 100000, 1}];
finalfunction = Interpolation[tabulate];
LogLogPlot[{finalfunction[y], y*finalfunction'[y]}, {y, 1, 1000000},  
  PlotRange -> {0.01, 10}, Frame -> True]

I used this same code on other similar solutions and had positive results. But for this specific solution I couldn't get it to work. Or it's seems that would take a very long time for it.

Could someone provide an alternative way to help Mathematica along?

Answer 1

If you break it up as below it will be considerably faster. The idea is to use symbolic integration on the inner integral (Why? Because we can, and it makes this faster).

In[263]:= 
i1[n_, t_] = 
 Integrate[
  Exp[-((z^2 + (200*n - z)^2)/(4*t))] + 
   Exp[-((z^2 + (100 + 200*n - z)^2)/200)], {z, -50, 50}, 
  Assumptions -> {t > 0, -4 <= n <= 4}]

Out[263]= (5*Sqrt[Pi]*(-Erf[10*n] + Erf[10*(1 + n)]))/
  E^(25*(1 + 2*n)^2) + 
   (Sqrt[Pi/2]*Sqrt[t]*(Erf[(25*Sqrt[2]*(1 - 2*n))/Sqrt[t]] + 
           Erf[(25*Sqrt[2]*(1 + 2*n))/Sqrt[t]]))/E^((5000*n^2)/t)

In[264]:= i2[t_] = Sum[i1[n, t], {n, -4, 4}];

In[265]:= 
pwd[y_] := (1/(4*Sqrt[Pi]))*
  NIntegrate[Exp[-1/(4*t)]/t^(3/2)*i2[t], {t, 0, y}]

The tabulation for the interpolation should complete in something approximating reasonable time. Well, within several hours, in any case, if counting by hundreds is any indication.

In[272]:= Timing[tabulate = Table[{y, pwd[y]}, {y, 0, 100000, 100}];]

Out[272]= {150.71, Null}