In the slow-motion region, ESR spectra cannot be expressed as a sum of simple Lorentzian lines. Studies of Freed and co-workers, on nitroxides in liquids gained information on the microscopic models of rotational dynamics, relying much on computer programs for simulation of ESR spectra based on the stochastic Liouville equation (SLE). However, application of Freed's method to copper system of biological interest has been for a long time precluded by lack of a full program able to simulate ESR spectra containing more than one hyperfine interaction. Direct extension of the Freed's approach in order to include superhyperfine interaction is not difficult from a theoretical point of view but the resulting algorithm is problematical because it leads to substantial increase in the dimensions of the matrix related to the spin-hamiltonian operator. In this paper preliminary results of a new program, written in C, which includes the superhyperfine interactions are presented. This preliminary version of the program does not take into account a restoring potential, so it can be used only in isotropic diffusion conditions. A comparison with an approximate method previously developed in our laboratory, based on a post-convolution approach, is discussed.