We derive cubic response functions in the Random Phase Approximation for calculations of non-linear frequency-dependent properties of molecules, and demonstrate an implementation of these functions using efficient computational algorithms. Illustrations are given by calculations of the frequency-dependent second hyperpolarizabilities, the frequency-dependent excited state polarizabilities and the three-photon transition amplitudes of the para-nitroaniline molecule, and by calculations of the frequency dependent magnetic second hyperpolarizabilities and their anisotropies for the first row hydrides isoelectronic with neon.