  We present a finite difference method to solve a new type of nonlocal
hydrodynamic equations that arise in the theory of spatially inhomogeneous
Bloch oscillations in semiconductor superlattices. The hydrodynamic equations
describe the evolution of the electron density, electric field and the complex
amplitude of the Bloch oscillations for the electron current density and the
mean energy density. These equations contain averages over the Bloch phase
which are integrals of the unknown electric field and are derived by singular
perturbation methods. Among the solutions of the hydrodynamic equations, at a
70 K lattice temperature, there are spatially inhomogeneous Bloch oscillations
coexisting with moving electric field domains and Gunn-type oscillations of the
current. At higher temperature (300 K) only Bloch oscillations remain. These
novel solutions are found for restitution coefficients in a narrow interval
below their critical values and disappear for larger values. We use an
efficient numerical method based on an implicit second-order finite difference
scheme for both the electric field equation (of drift-diffusion type) and the
parabolic equation for the complex amplitude. Double integrals appearing in the
nonlocal hydrodynamic equations are calculated by means of expansions in
modified Bessel functions. We use numerical simulations to ascertain the
convergence of the method. If the complex amplitude equation is solved using a
first order scheme for restitution coefficients near their critical values, a
spurious convection arises that annihilates the complex amplitude in the part
of the superlattice that is closer to the cathode. This numerical artifact
disappears if the space step is appropriately reduced or we use the
second-order numerical scheme.
