Generic Hamiltonian systems have a mixed phase space, where classically disjoint regions of regular and chaotic motion coexist. We present an iterative method to construct an integrable approximation H_reg, which resembles the regular dynamics of a given mixed system H and extends it into the chaotic region. The method is based on the construction of an integrable approximation in action representation which is then improved in phase space by iterative applications of canonical transformations. This method works for strongly perturbed systems and arbitrary degrees of freedom. We apply it to the standard map and the cosine billiard.