We report a numerical calculation of the two-photon absorption coefficient of electrons in a binding potential using the real-time real-space higher-order difference method. By introducing random vector averaging for the intermediate state, the task of evaluating the two-dimensional time integral is reduced to calculating two one-dimensional integrals. This allows the reduction of the computation load down to the same order as that for the linear response function. The relative advantage of the method compared to the straightforward multi-dimensional time integration is greater for the calculation of non-linear response functions of higher order at higher energy resolution. 78.20.Bh The measurement of two-photon absorption coefficient [1] yields different information from the single-photon absorption measurement, since the physical processes involved and selection rules are different. Because of this, there has always been a lot of interest in two-photon absorption of various molecules, crystals and solids [2][3][4][5][6][7][8][9]. To compare against the experimental data, one would like to have theorerical calculations based on some realistic model of the material. However, since the twophoton absorption is a typical non-linear optical process, its realistic modeling has always proved difficult for large complex systems [2][3][4][5][6].A powerful method that has come to be used widely for large quantum systems is the real-time real-space higherorder difference method [10,11], in which the real-space is represented by discrete mesh points, and the time development of a system is solved by numerically integrating the Schrödinger equation for discrete time steps. The energy levels and energy eigenstates are obtained by Fourier analyzing the numerical solution. The memory requirement scales linearly with the number of basis states N compared to N 2 for matrix diagonalization, and the method has proved effective in solving large quantum systems that cannot be solved by conventional methods [11,12]. So far, the large computation load has meant that the actual application of the method has been made primarily to the calculation of linear response functions of one-particle systems [10,11]. Nevertheless, the potential scale advantage of the method when applied to large systems invites the speculation that the development of the method will be essential in making the calculation of non-linear response functions of complex many-body quantum systems feasible.In this article, we report a trial application of the method to the calculation of the two-photon absorption coefficient of non-interacting electrons trapped by a binding potential and exposed to monochromatic light of frequency ω.The two-photon absorption coefficient is given by [1]where H is the unperturbed hamiltonian of the system, r is the electron coordinate operator, η is the frequency resolution, V is the volume of the system, the summation with respect to the initial state |E i must be taken over all states below the Fermi level E F , and the summation wit...