A numerical method for the solution of the time-dependent multigroup diffusion equations is presented. The method has the property that it is numerically unconditionally stable for all changes in reactor properties and all integration time-step sizes. The method assumes that the neutron flux and precursor concentration can be expressed as an exponential function over each time step. As a result of this assumption, and the factoring of the matrix form of the multigroup equations, it is shown that for the case of a constant step change in the properties of the system the asymptotic numerical eigensolution is proportional to the asymptotic eigensolution of the differential equations. An analysis of the truncation error associated with the method is also presented. Finally, a number of numerical experiments are presented which illustrate the accuracy, speed, and general utility of the method.