A steepest descent algorithm is used to obtain a least-squares approximation for threedimensional, toroidal stellarator flux surfaces represented by a discrete set of Poincare puncture points. A stream function \(p,0,0) is introduced as a renormalization parameter to improve the mode convergence properties of the double Fourier series for the inverse co-ordinates (R, Z) representing toroidally nested magnetic surfaces: R(p, 0,0) = 2 R mn (p) cos(m0 -n0) and Z(p, 0,0) = 2 Z mn (p) sin(m0 -n0), where p is a scaled radial flux co-ordinate and (R, 0, Z) are standard cylindrical co-ordinates. The variable 0 is the geometric toroidal angle, and 0 is a parametric co-ordinate representing a poloidal angle. The stream function X(p, 0,0) = 2 Xmn(P) sin(m0~n0) and the rotational transform profile t(p) are determined by solving a system of simultaneous linear equations obtained from the MHD equilibrium condition V X B-Vp = 0, together with the stellarator condition for zero net toroidal current on each flux surface,