And the same for plane = (6, 0.0). But I want to get a cross section when both the third variable (h_1) and the sixth variable (h_2) are zero at the same time.

Can I do it somehow? I tried to draw a plane through 2 six-dimensional points as in the example, but in this case the cross() function does not work.

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.

Something I didn’t understand from your first post was the goal of finding a Poincare section satisfying two conditions, h1(x) = 0 and h2(x) = 0. In general there’s no reason to expect that with only one independent variable (time t) you’ll be able to find points x(t) that satisfy two conditions simultaneously.

By definition, a Poincare section for the flow associated to a system of n first order differential equations is a (n-1)-dimensional surface (usually a hyperplane in R^n), that is transversal to the flow.
To choose a right Poincare section, not just a hyperplane of equation x_i =0 (or y_i=0 or zi_0) you should associate the 6d vector field,
V(x_1, x_2, y_1, y_2, z_1, z2), having as coordinates the right-hand sides of your system.

The general equation of a hyperplane in R^6 is Ax_1+B_x_2+Cy_1+Dy_2+Ez_1+Fz_2+G=0, and its normal is the vector, N, of coordinates [A, B, C, D, E, F]. If the dot product
of the vector field V and the normal, N, is nonzero
at any point of the hyperplane, <V(x_1, x_2, y_1, y_2, z_1, z_2), N> \neq 0, then that hyperplane is a right choice as a Poincare section.
Hence assigning particular values to A,B, C, D, E, F, you should find a Poincare section, if any.