Optical fields in metamaterial nanostructures can be separated into bright modes, whose dispersion is typically described by effective medium parameters, and dark fluctuating fields. Such combination of propagating and evanescent modes poses a serious numerical complication due to poorly conditioned systems of equations for the amplitudes of the modes. We propose a numerical scheme based on a transfer matrix approach, which resolves this issue for a parallel plate metal-dielectric metamaterial, and demonstrate its effectiveness.