We describe polarization models for Maxwell's equations which involve distributions of relaxation times and explain the difficulty in simulating these equations in the time domain for broadband pulsed forcing. We introduce a numerical approach based on Polynomial Chaos and discuss how to analyze this method.