You have a system of ODEs in the unknowns P(xi) and W(xi).  What is v that occurs in the second equation?  Saying that v is the velocity doesn't help.  It's the velocity of what?  Is it a constant?  If now, how is it related to P, W, and xi?

There have not been responses to your question after two days of posting.  That's because the code that you have posted is hard to read.  If you wrote that code in a Maple worksheet, you should post the worksheet itself.  Here is how.  In this web page go to the gray strip at the bottom of your question and hit "Reply".  That will bring up an editor where you may add a followup to your question.  In the editor's toolbar you will see a big fat green arrow.  Click on it and follow the instructions to upload your worksheet.

Furthermore, a part of your code reads:

# Define boundary conditions for the system
cond = 0, D(f1[0])(5) = 1]:
cond = 0, D(f2[0])(0) 1]: 
for j from 1 to N do
        cond = 0, D(f1[j])(0) = 0, D(f1[j])(5), D(f2[j])(0) = 0, D(f2[j])(5) = 0]:
end do:

That is broken beyond repair.  [Was it by chance written by Chat-GPT? ]  See if you can expand the comment by actually saying precisely what the boundary conditions are supposed to be.

@C_R It's possible to show that the function eta defined in my previous message is periodic if and only if the coefficient a is rational.  In particular, if a = k/m, where k and m are relatively prime integers (in other words, the fraction is reduced to the simplest form) then the period is 2*Pi*m.

In the animation, I picked a=1.3, that is a=13/10.  Therefore m=10, and eta is periodic of period 20*Pi.

To produce a non-periodic eta, pick an irrational number for a, such as sqrt(2) or Pi.

The reasoning is explained in periodicity.pdf.

@C_R The worksheet's angle psi measures the angle of the small circle's marker relative to the vertical.  As you have noted, the angle psi appears to be quite regular— its graph is something like that of the t + 4*sin(t) function — but that does not mean that the radius that connects the the circle's center to the marker will rotate in any regular fashion.  Here is a caricature of the situation.  Consider 

restart;
with(plots):
tmax := 17.86:
psi  := 1.3*(t + 4*sin(t));

The function psi generally grows with t indefinitely.  But if it represents a polar angle, we may mod out mutiples of 2*Pi from its value for practical purposes, which leaves us with

eta := psi - 2*Pi*trunc(psi/(2*Pi));

Here are what the graphs of psi and eta look like:

Note that eta is not periodic.  Here is what the radius in the direction eta(t) looks like in polar coordinates:

The erratic motion is similar to what we see in the rolling circles case.  I don't know much about chaos, so I can't tell whether what we are seeing is technically chaotic.

@C_R The rotation of the inner circle in the second example appears to be chaotic but I wouldn't know how to verify that.

@Earl The case where the outer circle is fixed is much simpler, but since it has only one degree of freedom, it's not a special case of the one treated here.  Doing it would be an instructive exercise.

@Kitonum That's a very clever solution.  I like it!

@vv Yes, you are right in that the Gateaux differential can be nonlinear in general.

In the problem at hand, however, the functional is Frechet differentiable if we assume q(w) is positive and bounded away from zero.  In that case the Gateaux and Frechet derivatives agree, and the differenial is linear.

@mayzal After making a few adjustments to your worksheet, I got pdsolve() to work and return a module. However, the evaluation of the module fails.  See the attached worksheet.

I have no further ideas on how to proceed.  Perhaps someone else may chime in.

question4-alt.mw

@mayzal One of your PDEs is

diff(X_dimensionless[1](t), t)  = ... P_dimensionless(z, t) ...

The right-hand side depends on z and therefore X_dimensionless[1] depends on z.  Shouldn't that be X_dimensionless[1](z,t)?

Aside: What is the point of the "_dimensionless" suffixes in your variable names?  Remove the horrible clutter by dropping those suffixes.  P(z,t) is as good a name for the dimensionless pressure as P_dimensionless(z,t), and is much easier to read.

@Kitonum The Wikipedia page Bernoulli quadrisection problem gives an example of an isosceles triangle which can be quadrisected in two different ways.

Kitonum, that's an excellent solution to a nontrivial problem.  Congratulations and vote up!

PS: In this post's first paragraph you write that your procedure applies to convex domains, but the examlpes show that it works with non-convex domains as well. Perhaps you meant simply connected?

@C_R I have no Maple Sim so I cannot be of help here, but let us know if you make further progress in what you are doing.

I don't have Maple Sim but I am guessing that what you are seeing is due to Maple Sim's attempt to simulate the dynamics of a massless disk.  Give the disk a little mass (something like 1e-6) and see if that helps.

Here is an interesting exercise: Calculate the vertical component of the force that the floor exerts on the disk. Is that force always possitive?

When that force becomes negative, the disk will lift off the floor!