The paper has considered an efficient method of pole shape optimization in multipole magnets to provide the given distribution of the magnetic flux density. The main idea of the proposed method is based representing the field characteristics in the form of a Fourier series and consequent compensation of undesired magnetic field harmonics by proper deformation of the pole tip surface. Fast convergence and observed robustness of the optimization procedure is explained by a strong correlation between the coefficients of the Fourier expansions for the magnetic flux density in the central part of the magnet and specific parameters used to describe the pole tip shape. The developed algorithm is integrated into computer software for simulating the magnetic field using the finite element technology developed for computing different types of electromagnetic systems.