@segfault You don't need to make a new Maple executable script.  As I wrote in my answer, define an alias in your unix shell.  If your shell is bash, for instance, then in your ~/.bashrc insert the line

alias maple='maple -j 65536'

This will take effect in subsequently opened terminals.  If your shell is something like tcsh, then in your ~/.tcshrc insert the line

alias maple  maple -j 65536

I am assuming that you launch Maple through the comand-line in a terminal.

@acer Thanks for replacing the page with a readable version.  I had noticed the issue with displaying vectors in prior posts and I had assumed that's due to a quirk of this website rather than related to the recent versions of Maple.

It's not your fault.  You calculations are correct.  I don't know why Maple returns FAIL. It looks like a bug to me.

If you replace the disk with a square, then Maple will return the correct result.

Aside: This is unrelated to the FAIL issue, but do you realize that the points that you are selecting in your code are not distributed uniformly on the disk?  They bunch up near the disk's center.

To get a uniform distribution on the disk, replace r1 and r2 with sqrt(r1) and sqrt(r2).

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?