  We present a numerical method for the computation of the conformal map from
unbounded multiply-connected domains onto lemniscatic domains. For $\ell$-times
connected domains the method requires solving $\ell$ boundary integral
equations with the Neumann kernel. This can be done in $O(\ell^2 n \log n)$
operations, where $n$ is the number of nodes in the discretization of each
boundary component of the multiply connected domain. As demonstrated by
numerical examples, the method works for domains with close-to-touching
boundaries, non-convex boundaries, piecewise smooth boundaries, and for domains
of high connectivity.
