In several previous publications we have outlined how a freeform optical surface can be tailored to achieve a desired irradiance on a reference surface. The shape of the surface is determined by numerically solving a differential equation. This approach has the advantage of being able to accommodate a huge amount of detail, equivalent to several thousand to millions of parameters for which classical optimisation techniques are not feasible. Laser beam shaping requires controlling the phase of the wave front as well. One surface is not sufficient, but two surfaces can be calculated by simultaneously solving a system of coupled differential equations. Justin L. Kreuzer patented this idea in 1965 (US Patent No. 3476463). Wassermann and Wolf outlined already in 1949 a procedure by which two aspheric surfaces are determined by solving differential equations such as to render an imaging optical system aplanatic. The work of Kreuzer as well as that of Wassermann and Wolf refers to rotational systems, where the cross section curve suffices to determine the surface and the equations are less demanding. We show how to extend the method of three-dimensional tailoring to simultaneously tailor two surfaces which are not necessarily rotationally symmetric, such as to achieve non-rotational irradiance distributions or off-axis devices. The desired irradiance translates into an equation which combines the first and the second derivatives of the surfaces.