A diagrammatic perturbation theory is presented which permits the calcn. of a mol. non-linear-optical response function for a particular time sequence of light pulses, or a mol. non-linear-optical susceptibility for a particular combination of frequencies of the input light fields. Mol. relaxation is treated in the Markov approxn. The method also accounts for incoherent population transfer between mol. levels as long as a hierarchical level ordering is possible. The diagrams are esp. useful for the selective generation of terms contg. a particular resonance of interest, or obeying the rotating wave approxn. The algorithm has been implemented in the symbolic algebra language REDUCE. These procedures allow the automatic generation of all diagrams evolving from a particular exptl. situation, and the automatic calcn. of the corresponding response function or susceptibility.