An easy way to do this is with polynomial interpolation. If your Poincare section is defined as h(x) = 0, integrate until h(x) changes sign. Then, for an nth-order integration scheme, make an nth-order polynomial interpolant of each the component of x as a function of h values, and interpolate to h=0.
If you use the Polynomial.jl package you can do this in a few lines of code. It’ll be allocating, and maybe there’ll be some cost to tracking h(x). Probably DifferentialEquations.jl has a slcik and efficient method.