Designing gradient coils for permanent magnets requires modeling the interaction between the currents and the ferromagnetic material. We have developed and tested an iterative design approach based on the boundary element method. The power dissipated by the gradient coil is treated as a constraint of the minimization problem. The technique is illustrated on a magnet with a rectangular prismatic cavity and two ferromagnetic plates adjacent to the coil. The procedure improved the linearity of the field from 10.5% to 3.4% and reduced the power dissipation to 63% of the power of the unoptimized coil. The optimized distribution of currents has been transformed to variable locations of coil windings without loss of gradient uniformity. The new approach significantly reduces computation time and memory storage and yields practical designs of gradient coils for permanent magnets.