Authors’ abstract: Based on the nonlinear time transformation method, in this paper we propose a recursive algorithm for arbitrary order approximation of heteroclinic orbits. This approach works fine for a wide class of systems that are perturbations of non-Hamiltonian integrable planar vector fields. Specifically, our method can provide an approximation up to any desired order, both for the locus where the heteroclinic bifurcation occurs in the parameter space and for the orbit in the phase space. We also give proofs that guarantee the existence and uniqueness of the solution. As illustrations, we apply it to the generic Hopf-zero and nonresonant double Hopf singularities which give rise to an intricate bifurcation scenario in the cases where a heteroclinic cycle appears in the corresponding unfolding of normal forms. The obtained high-order approximations are expressed in terms of symbolic nonzero normal form coefficients which improve the existing first-order approximations in the literature. They are also in excellent agreement with numerical continuations. Note that the derived formulas can be practically important in many control engineering applications in which certain constants of the problem may not be known in advance.