In the present paper we consider the efficient treatment of free boundary problems by shape optimization. We reformulate the free boundary problem as shape optimization problem. A second order shape calculus enables us to realize a Newton scheme to solve this problem. In particular, all evaluations of the underlying state function are required only on the boundary of the domain. We compute these data by boundary integral equations which are numerically solved by a fast wavelet Galerkin scheme. Numerical results prove that we succeeded in finding a fast and robust algorithm for solving the considered class of problems. Furthermore, the stability of the solutions is investigated by treating the second order sufficient optimality conditions of the underlying shape problem.