In this paper, one of the major shortcomings of the conventional numerical approaches is alleviated by introducing the probabilistic nature of molecular transitions into the framework of classical computational electrodynamics. The main aim is to develop a numerical method, which is capable of capturing the statistical attributes caused by the interactions between a group of spontaneous as well as stimulated emitters and the surrounding electromagnetic field. The electromagnetic field is governed by classical Maxwell's equations, while energy is absorbed from and emitted to the (surrounding) field according to the transitions occurring for the emitters, which are governed by time-dependent probability functions. These probabilities are principally consistent with quantum mechanics. In order to validate the proposed method, it is applied to three different test-cases; directionality of fluorescent emission in a corrugated single-hole gold nano-disk, spatial and temporal coherence of fluorescent emission in a hybrid photonic-plasmonic crystal, and stimulated emission of a core-shell SPASER. The results are shown to be closely comparable to the experimental results reported in the literature. arXiv:1907.09952v2 [physics.class-ph]