Why are my solution times not consistent with my results when using the solve command?

I'm trying to solve the heat equation for a 3-dimensional model of metal and high explosive. When the explosive reaches a certain temperature, it begins to generate heat exponentially. However, when I reach this point in my code, I receive the following errors:

Error using pde.TransientThermalResults (line 86)
The number of input solution times is not consistent with the results.
Error in pde.ThermalModel/solve (line 119)
    R = pde.TransientThermalResults(self,u,tlist);
Error in heat_diffusion_3D_laser (line 54)
S = solve(model,tlist);

If I select an end time for the solution ("Tmax") that terminates before this internal heating begins, the simulation performs well. The only explanation I can think of is that the internal heating increases the temperature for the solution so quickly that the times provided in "tlist" are insufficient for representing the gradients produced by solve. Code is below:

% casing properties (steel)
rho1 = 8050;   % density [kg/m^3]
k1 = 43;   % thermal conductivity [W/(m*K)]
cp1 = 502;   % heat capacity [J/(kg*K)]
L1 = 0.005;   % length of material [m]
% main charge properties (TNT)
rho2 = 7874;   % density [kg/m^3]
k2 = 80.4;   % thermal conductivity [W/(m*K)]
cp2 = 965;   % heat capacity [J/(kg*K)]
L2 = 0.05;   % length of material [m]
Q = 473;   % heat of decomposition [K]
Ea = 146e3;   % activation energy [J/mol]
Ze = 1.32e19;   % characteristic reaction rate (pre-exponential factor) [J/(kg*K*s)]
% misc properties
G = 8.314;   % ideal gas constant [J/(K*mol)]
T0 = 300;   % initial temperature [K]
Tmax = 9.82;   % engagement time [s]
% generate 3D model
modelr = 0.1;   % model radius [m]
gm = multicylinder(modelr,[L2 L1],'ZOffset',[0 L2]);
model = createpde('thermal','transient');
model.Geometry = gm;
pdegplot(model,'FaceLabel','on','FaceAlpha',0.5);
% assign cell properties
thermalProperties(model, 'Cell',1,'ThermalConductivity',k1,'MassDensity',rho1,'SpecificHeat',cp1);
thermalProperties(model, 'Cell',2,'ThermalConductivity',k2,'MassDensity',rho2,'SpecificHeat',cp2);
heatSourceValue = @(~,state) rho2*Q*Ze*exp(-Ea./(G*state.u));   % internal heat source due to HE heating
internalHeatSource(model,heatSourceValue,'Cell',2)
% assign boundary conditions
thermalBC(model,'Face',4,'HeatFlux',0);
thermalBC(model,'Face',3,'HeatFlux',0);
thermalBC(model,'Face',5,'HeatFlux',0);
thermalBC(model,'Face',2,'HeatFlux',@externalHeatFlux,'Vectorized','on');
% generate mesh
generateMesh(model,'Hmax',modelr/5);
figure;
pdeplot3D(model);
title('Model Mesh');
% initial conditions
thermalIC(model,T0);
% solve problem
tlist = linspace(0,Tmax,10);
S = solve(model,tlist);
% plot solution
figure
pdeplot3D(model,'ColorMapData',S.Temperature(:,end));
title('Final Temperature Distribution');

Also, the code for the external function:

function Qflux = externalHeatFlux(region,~)
P = 1450;   % beam power at target surface [W]
beamr = 0.025;   % beam radius [m]
A = pi*(beamr.^2);   % beam spot size [m^2]
I = P/A;   % irradiance at target surface [W/m^2]
R = 0.1;   % reflection coefficient at surface
[~,rho,~] = cart2pol(region.x,region.y,region.z);
Qflux = I*(1-R)*exp(-rho.^2/beamr.^2);
end