3-D graph of function with integral

Hi, I'm trying to graph a function of two variables containing an integral,

T(r,t)= (constants) * integral[ n^(-1.5) * exp(-r^2 / constants) ]

where the integral is evaluated for n, from bounds t-1 to t.

this equation models the temperature of a medium at a distance "r" and time "t" away from a point heat source. If I shift the timespan I'm graphing, that should not change the calculated temperature at a specific distance and time. But my surface plot shows some serious changes when I do... not cool. I've been trying to figure it out, but alas, I haven't found the key.

If anyone could shed some insight that would be fantastic!

Below is the code, with the relevant section after the note with all the asterisks. I just included the entire thing in case you want to see what I mean about the graphs: try one with

[r,t] = meshgrid(0:.000025:.001, 1.1:.5:21.1), and the other with

[r,t] = meshgrid(0:.000025:.001, 3.1:.5:23.1).

_________

%%Single Probe Point Heat Source
% T = f(r,t)

% establishing variables

tp = 1; % s
pwr= [0.0007184    0.0006090    0.0005146    0.0004355]; 
% ^ heat pulse power (per temp)

M = [
 273.15  999.84  4217.6  0.561  1.33E-07  -0.0181505  -0.0119493  -0.0183359  1938.97  2841.48  2163.05;
283.15  999.7  4192.1  0.58  1.38E-07  -0.0191278  -0.0125591  -0.0192895  1261.97  1829.71  1422.05;
293.15  998.21  4181.8  0.5984  1.43E-07  -0.0201162  -0.0131759  -0.0202541  821.53  1171.63  926.636;
303.15  995.65  4178.4  0.6154  1.48E-07  -0.0211047  -0.0137927  -0.0212187  541.691  760.925  608.064;   
];
%matrix M = [T rho Cp k a mR2 mR4 mR7 R2 R4 R7]

i=3; % run program for T=293.15 K (3rd matrix row)
T=M(i,1);
k=M(i,4);
a=M(i,5);
P=pwr(i); % W (depends on desired T though)

*********RELEVANT PART STARTS HERE*************

% WHY DOES CHANGING THE TIME BOUNDS CHANGE THE WHOLE PLOT??? (i.e. changes
% temp at the same time and radius)
[r,t] = meshgrid(0:.000025:.001, 1.1:.5:21.1);

% MUST MAKE SURE that "length" agrees with size(r) and size(t)
length=41
size(r)
size(t)

integ=[];

**********PROBLEM AREA HERE**************

% let's find the integral portion of the equation
for i=1:length % as r runs through 1:11
    for j=1:length % as t runs through 1:11
        rval=-(r(1,i)^2);
        fun=@(n) (n.^(-1.5)).*exp( rval *(n.^-1)/(4*a));
        integ(i,j)=quadgk(fun, t(j,1)-tp, t(j,1));
    end
end

T = (P/(8*pi^(1.5)*k*a^0.5)) * integ;
figure(25)
surf(r,t,T);
xlabel('Distance from Point Source (m)');
ylabel('Time Since Pulse (s)');
zlabel('Temperature Change (K)');