The implementation of quadratic response function in the four-component Hartree-Fock approximation was discussed. It was implemented based on the Dirac-Coulomb Hamiltonian. The second-order response function enabled the calculation of non-linear, time-dependent, corrections of the molecular polarization and magnetization. The first order response function was determined to evaluate the second-order response function.