This review provides a brief introduction to the physics of coupled exciton-plasmon systems, the theoretical description and experimental manifestation of such phenomena, followed by an account of the state-of-the-art methodology for the numerical simulations of such phenomena and supplemented by a number of FORTRAN codes, by which the interested reader can introduce himself/herself to the practice of such simulations. Applications to CW light scattering as well as transient response and relaxation are described. Particular attention is given to so-called strong coupling limit, where the hybrid exciton-plasmon nature of the system response is strongly expressed. While traditional descriptions of such phenomena usually rely on analysis of the electromagnetic response of inhomogeneous dielectric environments that individually support plasmon and exciton excitations, here we explore also the consequences of a more detailed description of the molecular environment in terms of its quantum density matrix (applied in a mean field approximation level). Such a description makes it possible to account for characteristics that cannot be described by the dielectric response model: the effects of dephasing on the molecular response on one hand, and nonlinear response on the other. It also highlights the still missing important ingredients in the numerical approach, in particular its limitation to a classical description of the radiation field and its reliance on a mean field description of the many-body molecular system. We end our review with an outlook to the near future, where these limitations will be addressed and new novel applications of the numerical approach will be pursued.