A volumetric finite-difference based method is presented in this paper for the design of three-dimensional (3D), arbitrarily structured gradient coils in magnetic resonance imaging. In the proposed method, the coil space is discretized with quasi-rectangular elements, and the current density of each element is expressed by a finite-difference based numerical approximation of stream functions. The magnetic flux density at target field points can be calculated by those stream function values at all grids of the coil space. The optimization problem is constructed and solved to determine the stream function and coil patterns. This proposed method has been tested on several designs that include a shielded, ultra-short cylindrical coil, a partially shielded biplanar coil, and an asymmetric head coil with 3D geometries. The numerical results show that the proposed method is straightforward to implement and is versatile and suitable for designing complex structured gradient coils with high electromagnetic performance.