Using the pseudo-spectral method for wave-packet evolution, we formulate and implement a new theoretical approach for studying the Bloch oscillation effect in a general number of dimensions. For this purpose, a novel recursive formula is obtained in the representation of Wannier functions for the infinitesimal time evolution of an electronic wave packet within the single-band approximation. This formula is valid for general position-dependent applied electric fields and arbitrary initial wave-packet shapes and crystal structures, and can also be generalized for time-dependent fields. For a tight-binding simple 'cubium' band, the infinitesimal time evolution can always be expressed in terms of Bessel functions of integer order; and for the 'empty-lattice model' band it can be expressed in terms of Fresnel functions for arbitrary electric fields also. As a further analytical application of our formalism, we derive an exact (finite-time) wave-packet evolution for the homogeneous and static electric field in the simple 'cubium' band. For this case, an estimate of the numerical error of the pseudo-spectral method is obtained for the mean position of the wave packet. As an illustrative example of the numerical implementation of our theoretical formalism, we present a simple one-dimensional numerical simulation for a Gaussian wave packet moving in a parabolic potential well. Finally, a general proof of the unitarity of the predicted time evolution is also presented, and different properties of our formalism pointed out.