Light incident on a layer of scattering material such as a piece of sugar or white paper forms a characteristic speckle pattern in transmission and reflection. The information hidden in the correlations of the speckle pattern with varying frequency, polarization and angle of the incident light can be exploited for applications such as biomedical imaging and high-resolution microscopy. Conventional computational models for multi-frequency optical response involve multiple solution runs of Maxwell’s equations with monochromatic sources. Exponential Krylov subspace time solvers are promising candidates for improving efficiency of such models, as single monochromatic solution can be reused for the other frequencies without performing full time-domain computations at each frequency. However, we show that the straightforward implementation appears to have serious limitations. We further propose alternative ways for efficient solution through Krylov subspace methods. Our methods are based on two different splittings of the unknown solution into different parts, each of which can be computed efficiently. Experiments demonstrate a significant gain in computation time with respect to the standard solvers.