@C_R,  well, I don’t know, I have Maple17, and it gives
only  
for both allvalues(solve(f, x, allsolutions)) and solve(f, x, allsolutions).

This is the equation of a fixed surface
f1 := x3-0.5e-1*exp(x1)/(0.5e-1+x1^4+x2^4);

The equation of a surface that rolls without slipping is shown in the picture.

And a simple torus with radiuses 1 and 0.5.  The picture shows the equation of the surface, rolling without slipping.

To clear my conscience, I added a spline connection to compensate for slight changes in the length of the middle link. This connection is implemented based on the polygon function.


OF_experimental_spline_connection_1_part_3.mw

@mmcdara  So it is, and so it is almost everywhere. I understand you perfectly.

@mmcdara 

 I regret that for many reasons I cannot maintain communication both on general topics and on professional issues. But I can enjoy watching such communication with your participation.
+ from me.

@Kitonum  
My attempts were on the forehead through the solution of a system of polynomial equations. I compiled several variants of systems of polynomial equations for the condition of the presence of an inscribed square, but they either could not be solved, or gave a very rough version of the solution.

@Rouben Rostamian  
We are talking about completely different dividers. I use the same normalization to move evenly along the curve, and I have the same length of 180. The message for vv was directly about the Draghilev method, and this is due to the possible absence of 0 in the denominator (the Cramer method of solving systems of linear equations is used there) when solving the ODE. And in the process of numerically solving the ODE, we can use normalization or not use it, that is, this has no direct relation to the Draghilev method.
Of course, if I understand you correctly.

@Axel Vogt 
Immediately I apologize for the quality of performance.
Part one. Here we show how a system of polynomial equations is obtained from the original system.
EQ_1.mw
Part two. Here we get a set of starting points for building solutions in Part 3.

In the third part, using the Draghilev method, we obtain one of the subsets of solutions, build graphs of its projections on our 3d, and selectively check the resulting solutions for the residual value. The solution period seems to be 2Pi.
EQ_2_and_3.mw

@Axel Vogt  Thank you for your interest. I have separate pieces of texts of the decision in rather amateurish execution. I will try to accompany them with comments.

@vv Draghilev's method is to get rid of the denominator in the ODE problem. This is the most fundamental. The Method works quite well without natural parameterization.
There are publications in Russian and English where the Method is applied. Some have good descriptions.

f1 := -8*cos(x1)+8*cos(x2)-8*cos(x3)-Pi*x4+4;
f2 := -2*cos(5*x1)+2*cos(5*x2)-2*cos(5*x3)+1; 
f3 := -2*cos(7*x1)+2*cos(7*x2)-2*cos(7*x3)+1;

This is a real practical challenge. It has an infinite number of solutions. It was solved by reducing to a polynomial system with subsequent accounting for periodicity. 
Draghilev's method was applied.

@Rouben Rostamian   Very nicely done and very well presented.

The intersection curve of the surfaces f1 and f2 is a connected and bounded set.
Moving along this curve, we track the points of its intersection with the surface f3.
These points are the solutions of our system. Number of intersection points 116.

@acer  Thanks, signchange=false is new to me.
(I just formally moved along the curve of the first equation and tracked the sign change of the second equation. There was a sign change at all 36 points.)

@Carl Love  
Didn't know about that, thanks. I found this site to help.
(I have a push-button phone and a PC. )