Abstract This research presents a new boundary integral equation method for approximating the conformal mapping w = ω(z) from bounded multiply connected circular domains of connectivity m + 1 in the z-plane onto the canonical domain obtained by removing m parallel slits from the strip − π2 < Im w < π2 in the w-plane. The method is based on a uniquely solvable boundary integral equation with the generalized Neumann kernel. The integral equation is solved using the Nystr¨om method with the trapezoidal rule. Several numerical examples are presented.