@Alfred_F 

I forgot to refer to "Niven I., Zuckerman H.S., Montgomery H.L. An introduction to the theory of numbers", page 183.

This problem is quite old. If you follow the sources in "Dickson", Volume 1, Chapter IX and the older sources that follow from it, then this and similar questions were already being dealt with around 1880. My attempt to use school-typical methods (factorization) to provide an understandable solution for amateurs was doomed to fail. Otherwise, mathematicians would have come up with it back then. The solution therefore requires more sophisticated number theory methods (Legendre). And one solution is shown in Carmichael's book "The Theory of Numbers", Mathematical Monographs No. 13, on pages 26 to 28, in principle using a slightly weaker, similar problem. The current problem here is given by Carmichael as exercise No. 9 with an asterisk (for a special level of difficulty).
Overall, your question is very interesting. It was a pleasure to study old sources. And in my view it is a gain to now see a solution formulated with the help of Maple's CAS.

@Ronan 

The task does not require that a solution be found in the Euclidean plane. There are surfaces in R^3 that are created by straight lines. And a quadrilateral with straight sides can lie on them. But that is then a very difficult task.

@Kitonum 

It would therefore be worthwhile to delve into the world of polygons. In a lattice plane (point lattice) Pick's theorem could then be of help.

@Kitonum 

Your program is instructive for me. I have now learned not only to be able to read command structures, but also to understand them.

The theoretical/topological background is well known. The quoted sentence only proves the existence of area division. You have gone a step further and constructed a solution. This raises the question: Are there further divisions for the same Jordan curve? At what point in your program would something have to be changed?

It is worth trying to substitute: ln(y(x))=z(x). The left-hand side of the differential equation becomes z´(x) and for x>1 the right-hand side is continuous. Solvable problems of a similar structure can be found, for example, in the terrifying collections:
- EqWorld
- Kamke differential equations. Mathematica 12.3 and Maple 2021, Nasser M. Abbasi
- https://dlmf.nist.gov

@Ronan,  @vv

I found this task while cleaning up one of my hard drives. A long time ago, with some help, a solution was found in Mathcad 14 code. This solution is very long, but it makes the way of thinking on the way to the solution clear even for interested amateurs. That was my goal at the time. Unfortunately, I can only do simple tasks in Maple so far. But the solution presented here helps again. The brevity in dealing with ODE and graphics is impressive.

@mmcdara 

So the small orthogonal lines mean that all sides are the same length. That was new to me.

According to Picard Lindelöf's theorem, the interval (approximately lying at) -0.5<x<1 must be considered separately. Here there are discontinuities on the right-hand side of the ODE. Continuity and the existence of the solution are only guaranteed outside of it.

What do the short orthogonal crossbars on the polygon sides mean? I don't know this symbolism. Are the polygon sides all the same length? In this case, the task would be easy to solve. Otherwise, at least ABC forms a kinematic chain. Further specifications would be required for a clear solution.

My first answer needs to be corrected. I didn't notice that the ODE was presented in an abbreviated form in the posting. I therefore missed the square of the first derivative.

...there under number 2.422 on page 504 you will find the solution: y=A*cot(x)+B*(1-x*cot(x)). Substitution according to Riccati y´´=y´*u could help.

@vv 

... and I wish you a peaceful Christmas.

@vv 

I know this element-wise proof. I had hoped to find hints for a different solution in the "Group Theory" package of the help text. To do this, it would be necessary to use the symbolic calculation rules for regular matrices to solve the equation A*B = B*A step by step for B without element-wise calculation. In this, B is the matrix we are looking for. One approach to this would be B-A^(-1)*B*A=0 and

det(B - A^(-1)*B*A)=0. With x for the eigenvalues ​​and regularity of A^(-1)*B*A, det(E*x - A^(-1)*B*A) = 0 would also be possible. From this we can guess the well-known result.

@Kitonum 

The solution is the period of the fraction 1/7.