A general density matrix treatment is presented for the localized dynamics of a molecular system strongly coupled to a medium, and leading to dissipation and fluctuation phenomena. A self-consistent field description allows for the response of the medium when this is driven away from equilibrium. The dynamics of the primary region is described with a dissipative potential within an atomic model, whereas the medium is described in a statistical manner in terms of dissipative rates. This treatment is applied to the femtosecond photodesorption of CO from the Cu(001) surface, in a model which includes several vibrational modes of the adsorbate, and the nonlinear response of the substrate metal treated with a modified optical Bloch equation. A computational procedure based on the split operator propagator and fast Fourier transform is applied to several studies with an increasing number of adsorbate degrees of freedom. It introduces a combination of discretization on a grid and expansions in basis sets of vibrational adsorbate states, to obtain results on yields of desorbed CO versus the fluence of the light pulse, in very good agreement with experimental results.