The aim of this paper is to describe a Matlab toolbox, called µ-diff, for modeling and numerically solving two-dimensional complex multiple scattering by a large collection of circular cylinders. The approximation methods in -diff are based on the Fourier series expansions of the four basic integral operators arising in scattering theory. Based on these expressions, an efficient spectrally accurate finite-dimensional solution of multiple scattering problems can be simply obtained for complex media even when many scatterers are considered as well as large frequencies. The solution of the global linear system to solve can use either direct solvers or preconditioned iterative Krylov subspace solvers for block Toeplitz matrices. Based on this approach, this paper explains how the code is built and organized. Some complete numerical examples of applications (direct and inverse scattering) are provided to show that µ-diff is a flexible, efficient and robust toolbox for solving some complex multiple scattering problems.