var('t')
PThStats = plot(thMean(t), (t, tMin, tMax), color = 'black') + \
plot(thMean(t) + thStD(t), (t, tMin, tMax), color = 'black', linestyle = ':') + \
plot(thMean(t) - thStD(t), (t, tMin, tMax), color = 'black', linestyle = ':')
show(PThStats, axes_labels = ('$t$', '$N$'))