@C_R  

The Isolate procedure finds all real solutions to polynomial systems of equations. In our case, the system of polynomial equations is the original system + f8 and f9. f8 sets the trajectory, f9 helps find the starting point for the Drag(h)ilev method. The original system, that is, without f8 and f9, has an infinite set of solutions, which, in turn, consists of an infinite number of infinite subsets.In the figure, we can only attempt to depict a single subset of solutions. In this case, it satisfies both the tangency condition and the condition that all spheres are inscribed between the surfaces, i.e., for the case where i=6.

Now try replacing the corresponding lines in the text with the following:

R := 2; r := 1; h := 0.2e-1; smin := -2.; smax := 8.;
------------------------------
T := Isolate([f1, f2, f3, f4, f5, f6, f7, f8, f9], [x1, x2, x3, x4, x5, x6, x7, x8, x9]): 
j := nops(T): i := 4;
-------------------------------
U2 := sphere([.5, 0, 0], r, transparency =0.6, color = [yellow], style = surface):

And you'll see that the "inscribed" sphere will touch the cylinder and absorb the stationary sphere. The tangency condition will be met, but the inscription condition will not. That is, this is a different subset of solutions—not ours. Approximately the same thing can happen with other values ​​of i. This can easily be determined manually using  simple graphs, since i istypically finite and small. 
In the case of transcendental equations, it's a bit more complicated, but essentially unimportant. For example, you can use the fsolve procedure with intervals to limit the search space.

@C_R 
Thank you for such attentive attention to my message. I didn't expect that. Well, there's nothing I can do about it for now. I really hope it doesn't interfere with understanding of the algorithm. In any case, I'm always happy to answer any questions.

@C_R  

The drive can be implemented in various ways, but from what I've heard and seen, lever mechanisms typically rely on uniform input rotation. Incidentally, I can't find one right now, but there is a video showing uniform rotation in the case of a Schatz mechanism. Apparently, this can be achieved through a regular gearbox. 
If you ever need this, the corresponding normalization for solving ODEs in the Draghilev method allows you to do so. This is included in the text.

@C_R 
Hi. Well done! I looked, but not everything. We have significant internet restrictions at the moment. I even understood some of what you wrote.
But I couldn't figure out for sure whether the input link rotates uniformly. It's not obvious to the naked eye, but your word will suffice. (Here, - in the topic, in all messages the left struts rotates uniformly.)

@vv 
Yes, the last option is the solution.

@C_R 
"Does this mean that the trajectory is given? By that, do you mean the trajectory of the contact point?"
Yes, of course. These are just graphs that appear in turn, reminding  of kinematics.
A small curved torus (or paraboloid) is tied with a set of points on its surface (intersection line), this set of points rolls along another set of points (trajectorie) on a stationary surface. In each frame, the equation of the moving surface is transformed based on the position of the points, and we get another graph and, accordingly, a frame.

@Rouben Rostamian  
You explained much better than me, and in a few words, how the algorithm works.
And special thanks to you for the  "Драгилев".

@Rouben Rostamian  
In several of my old posts there was a proposal to apply the Method to solving the inverse kinematics of manipulators (if you want, I will try to find all these messages). There we proceed as follows: we impose additional connections on the links of the manipulator. We transform it "mathematically" into a spatial lever mechanism with one degree of freedom. After which the inverse problem is solved automatically. And we choose the type of constraints from the simplest ones, so that the working point can easily move along its trajectory.
Here we use a similar technique for rolling a figure on a surface. For example, we take a small torus, intersect it with a plane. In the same way, we obtain a trajectory on a stationary surface, intersecting it with some surface. We make the step along both intersection curves the same there and there. We tie the equation of the small torus to the points of the section line and roll the curve along the curve. But we do not show these curves, we only show the graphs of our surfaces. The plane of the points of the section curve on the small torus is directed approximately along the fixed curve in the osculating plane of this fixed curve and is additionally corrected by the normal to the fixed surface at each point. We can say that we leave it 1 degree of freedom, which is what we see in the figure.
(I keep using the letter "h" in his last name, but now I don't know if it's right or not.)

@C_R  
Did you like it? Thank you, I thought a lot on the algorithm. In fact, nothing rolls in my algorithm. The rolling surface is in its initial state before each new frame. Then we look at where its current location on the trajectory should be, based on the very first point of the trajectory,
we apply Euler angles with an offset, and after memorizing the corresponding graph, we send the rolling surface to its initial state. And from which side to roll, we determine using the normal to the stationary surface.
Honestly, I have not yet encountered such 3D animations with surfaces, and it would be interesting to see rolling surfaces without slipping performed by MapleSim.
In any case, these are just toys. The only thing I will draw your attention to is that the Dragнilev method is also used here.☺

And just a little bit more

Adding images on this topic

@acer 
Are you talking about G1[iii]? I couldn't think of anything else. After all, each frame produces a new equation, and it needs to be displayed on the graph. And even then, the animation is shown after 2 frames due to resource limitations.

@Rouben Rostamian  No, it's much simpler. This is computational geometry. A set of points is rigidly attached to a curved (small) torus, and a trajectory is laid out on the large torus, along which the small torus supposedly rolls. At each point of the trajectory, calculations are made so that the small torus imitates rolling without slipping. This is done using Euler angles: we recalculate the coordinates of the attached points, and then, based on this, we transform the equation of the small torus and draw another graph of motion.
I'll never be able to handle physics equations.

@Rouben Rostamian  
I thought  if we move along the graph line itself, then the proximity between the roots has a weaker effect on the situation with their possible omission. That is, we are talking about the values ​​y(s) and x(s), where s is a natural parameter (arc length). And then we can, for example, use a cycle "to count" harmonics based on their average length. And each time start counting from the last point on the curve. Then we should equally successfully pass "empty" ones, where y does not change sign, and places with a cluster of roots. In this case, most likely, it is necessary to additionally clarify each last point on the line at each step of the supposed cycle. 

(edited)
No, it looks like it won't be more reliable than your proposal. To move along the curve, an unacceptable number of significant digits would be required. "My" method is already very poorly controlled beyond 10^5

@Rouben Rostamian  

I don't need all the solutions of the cluster. I want to make sure that at a significant distance, for example, at a distance of 10^5 - 10^6 from the origin, not a single root out of about 50 roots in a row is guaranteed to be skipped. I'm just comparing with the capabilities of that very method of continuation by parameter. Here, I recently picked up an example. I'm already embarrassed that I attracted so much attention to myself, but it is difficult for me to satisfy my curiosity on my own due to my knowledge of the package.