@hamideh 

The whole animation is a solution on a tube; the solution can be continued in a practical sense for the entire tube. This is not one or three points, it is an infinite number of solutions. Look carefully at the text. This method can find all solutions on connected sets (for example, on all our tubes separately). The method only needs to hook on one point of such a set. In order to catch on, you need all sorts of techniques, including auxiliary equations.
I advise you to familiarize yourself with the idea of the method. You can do this here.
https://www.maplesoft.com/applications/view.aspx?SID=149514
A lot of information in Russian. There is also a search on the forum. To better master the Draghilev method and help others with this, you can create a separate topic in the questions section.
The algebraic solution is hardly better than the solution by Draghilev's method, but in this case, it seems to me, is much more convenient. By the way, it was found with help by Draghilev's method.
By controlling three parameters: first the radius of the cylinder, then two angles of inclination, you can very accurately select the geometric appearance of the first tube. Then two more parameters remain: the period along the first coordinate and the period along the second coordinate. The quality of the selection of all parameters can be checked by the accuracy of the difference between the original equation and the new one.

@hamideh 
The solution can be represented by an infinite number of inclined cylinders with a radius of Pi / 20.
Perhaps for this it is only necessary to choose the angles of the rotation matrix a little more precisely.
The approximate periods of the coordinates are, respectively, (+ -2 * k * Pi / (7.7)) and (+ - k * Pi / (6.7))).

restart; with(plots): 
r := (1/20)*Pi; 
a := 0; 
b := evalf(Pi/(5.1));
c := evalf(-Pi/(4.)); 
m[0, 0] := cos(a)*cos(b); 
m[0, 1] := (-cos(a)*sin(b))*(-sin(c))+sin(a)*cos(c); 
m[0, 2] := -cos(a)*sin(b)*cos(c)+sin(a)*sin(c);
m[1, 0] := -sin(a)*cos(b); 
m[1, 1] := -((-sin(a))*(-sin(b)))*sin(c)+cos(a)*cos(c); 
m[1, 2] := ((-sin(a))*(-sin(b)))*cos(c)+cos(a)*sin(c); 
m[2, 0] := sin(b);
m[2, 1] := -cos(b)*sin(c); 
m[2, 2] := cos(b)*cos(c); 
xn1 := x1*m[0, 0]+x2*m[0, 1]+x3*m[0, 2]; 
xn2 := x1*m[1, 0]+x2*m[1, 1]+x3*m[1, 2]; 
xn3 := x1*m[2, 0]+x2*m[2, 1]+x3*m[2, 2];
 f1 := -r^2+xn1^2+xn2^2;
 f2 := cos(2*Pi*(x1+x2-2*x3))+cos(2*Pi*(x2+x3-2*x1))+cos(2*Pi*(x3+x1-2*x2));
 f3 := (xn1-2*Pi/(7.7))^2+(xn2+Pi/(6.7))^2-r^2; 
Sf := implicitplot3d([f1, f2, f3], x1 = -.75 .. .75, x2 = -.75 .. .75, x3 = -.75 .. .75, color = [RGB(.5, .1, .21), blue, green], numpoints = 5000, transparency = .5, style = surface); 
display(Sf, axes = normal)

It turns out that in the end we can work with an algebraic equation, the equation of a cylinder.

 

@hamideh 
In the first part of the program, we build the curve of intersection of one of the "tubes", for example, with a plane. In the second part of the program, we build curves that intersect our "tube" perpendicular to the first curve at each its point. 
At the end of the work of each part of the program, the result is displayed in the form of a graph. The last graph is animated.
The set of arrays Lg [j] (j = 1..3) is a numerical solution of the selected subset of all  (x, y, z) on a particular "tube".
The distance between the centers of the "tubes" seems to be about  3*(1.399 / (2*Pi)).
The number of solution points and their accuracy can be controlled using parameters in the numerical solution of autonomous ODEs.
If you need to improve the accuracy of the solution, you can use, for example, Newton's method at each point. But it will require additional counting time.
I ask you to consider this text of the program solely as an idea of the algorithm, but not as a program written in Maple. 
F(X_Y_Z)=0_hamideh.mw

@hamideh 
The graph is not a solution; rather, it shows what the solution looks like. I will find the time and try to show the application of the Draghilev method for your task. This should be one of the "tubes" in numerical form. Then from this set of coordinates you yourself will be able to choose the ones you need taking into account of the period for all the others.

@Kitonum  Гениально, Юрий Николаевич. 

@vv  the numerical solution of an ODE cannot claim to be accurate when finding solutions to nonlinear equations. Each such task must be approached individually. 

@vv  @mmcdara 

This is a collaboration of the Draghilev method with NextZero. I ask you not to pay attention to the quality of the text.
A new variable v is introduced into the original equation; at its zero value, we get a solution to the original equation.
x1*(1+1.000100000*sin(x1^2)) - 2.129280519*v = 0;
In the program, we track when this variable takes a value of 0 and fix the solution to the original equation.
The starting point of the original variable in this case is 1.1, and the initial value of v is always 1 (this is like in the method of continuation of the solution by parameter).
For example, for the integration interval from -50 to 0 and for the starting point 1.1, 3614 solutions are found in the positive direction (printed in reverse order due to the integration interval).
Draghilev_NextZero.mw

@mmcdara  Come on, it's okay. The main thing is that we figured out the reason.

@vv  Thank you for your so kind feedback. 

@mmcdara  Look where the function changes sign

restart; 
x := (1/2)*(4.155071365+4.158473414); f := x*(1+1.0001*sin(x^2));
 x := 4.155071365; f := x*(1+1.0001*sin(x^2));
 x := 4.158473414; f := x*(1+1.0001*sin(x^2));

All solutions are additionally verified by Draghilev's method. If interested, this method didn't miss any solutions.

@Carl Love  Think this coefficient is not directly related to Draghilev's method. I don’t know how to get a closed form solution for this autonomous differential equation system (and Maple does not show the progress of the solution). The expression will probably look better if there is no such factor for s. But I don't know how to do it. I only have experience with parameterizing the arc length of a curve for a numerical solution.

@Carl Love  I just have no words. Thank you. You and vv show me what Maple is.  

Here 

It seems that Draghilev's method works more "neatly". This can be seen in the example of a line segment with ends
[0.5, 1, 1.5], [-1, 0, 1]. If we weaken the attention to the solution based on Gröbner bases and do not additionally divide the segment into parts, then formally possible, but in reality hardly admissible discreteness may arise: 
This, as we can see, does not happen when using Draghilev's method:

@Ronan   
An example of exactly the same movement as in the first example  PLAT 1. To do this, replace the third part of the program with this text:
STEWART_PLATFORM.mw
Platform design options are endless.