The Levi-Civita wavefront coupling problem consists of, given two prescribed wavefronts, obtaining a refractive or reflective surface coupling them. We propose a robust numerical method to solve Levi-Civita's problem, whose rationale is to consider that Levi-Civita's solutions are level surfaces of the point characteristic function established between points of incoming and outgoing wavefronts. The method obtains both surface data points and their normals, enabling a more robust surface reconstruction. We carry out a detailed error analysis of our method by means of comparing the surface data estimation with nominal surfaces obtained in reference tests offering analytical solutions to Levi-Civita's problem. The method offers, in computer simulations, highly accurate results with moderate computational cost.