A real-space method is developed to calculate molecular hyperpolarizabilities in the time-dependent density functional theory. The method is based on the response function formalism which was developed by Senatore and Subbaswamy for the third harmonic generation of rare-gas atoms [Phys. Rev. A 35, 2440 (1987)]. The response equations are discretized in real space employing a uniform grid representation in the three-dimensional Cartesian coordinate, and are solved with iterative methods such as conjugate-gradient and conjugate-residual methods. The method works efficiently for both small and large molecules, and for any nonlinear optical processes up to third order. The spatial convergence of the calculation can be examined with two intuitive parameters, the grid spacing and the spatial box size. Applications of our method are presented for rare-gas atoms and molecules, N-2, H2O, C2H4, C6H6, and C-60. Our results agree well with other calculations employing basis functions except for a slight deviation in a large molecule, C-60.