It is still very challenging to determine a freeform lens for converting a given input beam into a prescribed output beam where not only the irradiance distribution but also the phase distribution hardly can be expressed analytically. Difficulties arise because the ray mapping from the input beam to the output beam is not only intertwined with the required double freeform surfaces but also intertwined with the output phase distribution, whose gradient represents the directions of the output rays. Direct determination of such a problem is very difficult. Here, we develop a special iterative wavefront tailoring (IWT) method to tackle this problem. In a certain iteration, the current calculation data of the double freeform surfaces and the output phase gradient are used to update the coefficients of a Monge–Ampère equation describing an intermediate wavefront next to the entrance freeform surface. The solution to the wavefront equation could lead to an improved ray mapping to be used to update the corresponding phase gradient data and reconstruct the double freeform surfaces. In a demonstrative example that deviates much from the paraxial or small-angle approximation, the new IWT method can generate a freeform lens that performs much better than that designed by a conventional ray mapping method for producing two irradiance distributions in the forms of numerals “1” and “2” on two successive targets, respectively.