We present a numerical study on the spatial distribution of fluorescence and photobleaching occurring in samples subject to multi-photon excitation. We developed a simulation model and implemented a simulator program. Its quantitative predictions can help to find the optimal operating parameters (such as laser power, pulse length, pulse repetition rate) of the two-photon microscope to reach higher image quality, to reduce undesired photobleaching, and to pave the way for optimized photoswitching-based super-resolution imaging. Conversely, the simulator might also be useful when photodynamic parameters are searched for. Furthermore, such simulations can promote the evaluation of the results of other fluorescence-based techniques [e.g. fluorescence recovery after photobleaching (FRAP) measurements]. The photodynamic model of the fluorophore contains a ground state, an excited state, a triplet state, and several photobleached states; the state transitions are characterized by absorption cross sections and lifetimes. The sample is modeled as a fluorophore solution divided into cubic cells among which diffusion takes place. The illumination is simulated as a focused laser pulse train described by a pulsed Gaussian beam. As a demonstration of the capabilities of the simulator, an example is presented that reveals the spatial distribution of photon emission in the sample investigated by a two-photon microscope in the case of different laser and photobleaching parameters, assuming one-photon absorption induced photobleaching. The simulation demonstrates quantitatively how photobleaching affects the spatial distribution of fluorescence and the resolution of the microscope.