So, what's your theorem, or the featured item in your paper? Is it that the maximal gap for odd primes p and q is min(p, q) - 1?

@Alfred_F I don't know if you consider the following to be a "trick", but it certainly still needs to be done to complete @mmcdara 's solution: Generate the 2nd initial condition. It can be done like this:

IC_2:= eval[recurse](convert(eq, D), {x= 0, IC});
                                        3
                      IC_2 := D(y)(0) = -
                                        2

@Ronan The Beta Forum is a forum for beta testers of the next-to-be-released version of Maple. See https://beta.maplesoft.com/welcome/ and https://beta.maplesoft.com/about/

One needs to apply and be approved to get access. Details are on the above web pages.

While cleaning up that worksheet for posting, I accidentally deleted a line of code. (I need to acquire the habit of always doing "Execute entire worksheet" before posting.) The corrected worksheet is now above.

@WA573 I think that the command that you want is VariationalCalculus:-EulerLagrange.

@C_R You might not realize that you are a Moderator, and thus you can read the deleted content by attempting to Edit it. If you do that, you will see that all 12 of those deleted Answers are pure spam: advertising for something related to the game Minecraft. The vast majority of moderator-deleted content is spam, as you can verify.

I don't think most readers here realize the huge amount of spam that I and other Moderators delete every day. By "spam" I mean strictly unsolicited advertising for products and services that have absolutely nothing to do with Maple or math. I never use the word to refer to any other type of low-quality, inappropriate, or repetitive writing, no matter how much I think it should be deleted.

@dharr It's a shame that some of Maple's integer-arithmetic commands (such as isqrt and iroot) return the closest value rather than the floor. I've been computer programming for 46 years. When doing exact integer arithmetic, I almost always want the floor, I occasionally want the ceil, and I've absolutely never had an application where I wanted the closest value. So, every time I use iroot or isqrt, I need to adjust the result.

Can anyone recall ever wanting the closest exact-integer value?

Is an annulus allowed as one (or more) of the parts? Or are only topological homeomorphs of a disk allowed? If annuli are allowed, is their circumference defined to be the sum of the circumferences of their inner and outer circles? 

@dharr I just posted a corrected Answer. Although I could've used shape= hermitian, this didn't seem to make any noticeable improvement in the Eigenvectors performance.

@dharr Egad, you're right. Sorry, about that. I'm going to try with shape= hermitian instead.

@Susana30 I added numerous "bells & whistles" to the plot to include features from your hand-drawn one.

plots:-display(
    plots:-polarplot~(
        [r1, r2, [r, theta__1, r= 0..6]], theta=~ [Pi/2..Pi, (0..Pi)$2], 
        color=~ [red, blue, green], thickness=~ [2,2,3],
        legend=~ [r[1]=r1, r[2]=r2, theta[1]=theta__1], 
        axis[2]= [tickmarks= 13], axesfont= [times, bold, 12],
        labels= [`\nr`, theta], labelfont= [times, bold, 18],
        coordinateview= [default, 0..Pi]
    ),
    plots:-polygonplot~(
        [Circle__arc, Cardioid__arc], transparency= .5,
        color=~ [pink, "LightBlue"], legend=~ ["Area 1", "Area 2"]
    ),
    legendstyle= [location= left, font= [times, 14]], size= 500*[2,1]
);

@Muhammad Usman Is it possible that in the example that seemed to work, the graphs were symmetric with respect to the main diagonal, the line y=x? If so, that would explain why it worked.

@Carl Love Here are two more algorithms for the sort. I've included a threadsafe version of degree, which makes each procedure below threadsafe in its entirety. If you don't need it to be threadsafe, you can replace degree__TS with degree, which'll probably be faster. For both algorithms, Nterms is the same as above.

The first new algorithm partitions the list by the signifying variable, then sorts each block by degree and number of terms, then the sorted blocks are catenated. If there are any entries that contain none of the variables. they are sorted by number of terms and placed at the end:

#Thread-safe version of degree, optimized to only be used on expanded 
#polynomials that contain v:
degree__TS:= (p,v)-> max(1, op~(2, indets(p, identical(v)^anything))[])
:
PolySort2:= (F::list, V::list(name))->  local H, N:= F, v;
    [
        (
            for v in V do
                (H,N):= selectremove(has, N, v);
                sort(H, 'key'= [rcurry(degree__TS, v), Nterms])[]
            od
        ),
        sort(N, 'key'= Nterms)[]
    ]
:

The second algorithm also partitions by the signifying variable, then each block is sorted into bins by degree of that variable, then each bin is sorted by the number of terms. Finally, it's all catenated. Entries that contain none of the variables are treated as above: 

#Variant of ListTools:-Classify optimized to this purpose. Returns a 
#list of lists instead of a table of sets.
Classify:= proc(f::appliable, L::list) local e, T:= table();
    for e in L do T[f(e, _rest)][e]:= () od;
    [entries](indices~(T, 'nolist'), 'indexorder')
end proc
:
PolySort3:= (F::list(algebraic), V::list(name))-> local H, N:= F, v;
    [
       (
            for v in V do
                (H,N):= selectremove(has, N, v);
                (op@sort)~(Classify(degree__TS, H, v), 'key'= Nterms)[]
            od
        ),
        sort(N, 'key'= Nterms)[]
    ]
:
PolySort3(F, [x4,x3,x2,x1]);

If you have large lists to be sorted, I'd be curious which algorithm is fastest.

@emendes 

1) Nterms is not an invariant:

The number of terms in an algebraic expression is a "structural" or "syntactic" property of the expression, not an algebraic or mathematical invariant. Your simplify(F1) is partially factored. There's no practical way of counting terms without expanding the polynomial (due to the possibility of cancellation of terms). That's why I said earlier that I assumed that all the polynomials were expanded. The number of terms of an expanded polynomial is an invariant. If you need to work with unexpanded polynomials, it'd be easy to modify the code below to expand each polynomial before assigning it a key. The expanded form needn't be stored. However, the time for expand could easily be the majority of the time.

2) Passing the variables to the Key:

Yes, it can be done. Here's some code:

Nterms:= p-> `if`(p::`+`, nops(p), `if`(p=0, 0, 1))
:
PolySort:= (F::list, V::list(name))->
local S:= {F[]}, B:= 1+max(map(op@degree~, S, {V[]})[], Nterms~(S)[]);
    sort(
        F,
        'key'= proc(p) local v, d, r:= 1;
            for v in V do d:= degree(p,v); r:= B*r+d until d<>0; 
            B*r + Nterms(p)       
        end proc
    )
:
PolySort(F, [x4,x3,x2,x1]);

This code is not threadsafe because degree is not threadsafe.

@Carl Love The list-form keys above help to understand how keys and lexicographic sorting of lists works. But we can easily make the key an integer, which is likely more efficient. Just consider the list elements to be digits of a base-B number where some sufficiently large number is substituted for infinity, and the base becomes infinity+1. To be "sufficiently large", the number just needs to be larger than any other digit that might appear. So, in this case, it needs to be larger than both the largest single-variable degree and the largest number of terms. The procedure below uses base B:= 100, mostly because that makes it easy to see how the keys are obtained. For this example, B:= 7 would work just as well.

Nterms:= p-> `if`(p::`+`, nops(p), `if`(p=0, 0, 1))
:
Key:= proc(p) local v, d, B:= 100, r:= 1;
    for v in (x4,x3,x2,x1) do d:= degree(p,v); r:= B*r+d until d<>0; 
    B*r + Nterms(p)       
end proc
:
sort(F, key= Key);

Key~(F);
             [10206, 10103, 10104, 1000206, 1000205, 100000205]