The construction of C ² Pythagorean-hodograph (PH) quintic spline curves that interpolate a sequence of points p ₀,...,p _N and satisfy prescribed end conditions incurs a “tridiagonal” system of N quadratic equations in N complex unknowns. Albrecht and Farouki [1] invoke the homotopy method to compute all 2 ^N+k solutions to this system, among which there is a unique “good” PH spline that is free of undesired loops and extreme curvature variations (k∈{−1,0,+1} depends on the adopted end conditions). However, the homotopy method becomes prohibitively expensive when N≳10, and efficient methods to construct the “good” spline only are desirable. The use of iterative solution methods is described herein, with starting approximations derived from “ordinary” C ² cubic splines. The system Jacobian satisfies a global Lipschitz condition in C ^N , yielding a simple closed-form expression of the Kantorovich condition for convergence of Newton–Raphson iterations, that can be evaluated with O(N ²) cost. These methods are also generalized to the case of non-uniform knots.