Search for "color bar" in the MaplePrimes search. Many solutions have been posted over the years. There is currently no formal mechanism for these color bars.

@Hullzie16 You can include your inequalities and also use the options parametric, real:

solve({f=0, discrim(f,y)>0, m<0, z>0}, {y}, parametric, real)

This reduces the number of cases in the piecewise. In this instance, they're reduced to just 1 case (not counting the "otherwise", which doesn't satisfy the inequalities anyway). 

@jud Okay, given the different way that the labels are placed in the graph in Maple 2021 (raw integers rather than strings), the following will work for you:

restart:
Gp:= GroupTheory:
S4:= Gp:-SymmetricGroup(4):
L:= Gp:-GroupOrder~(Gp:-SubgroupLattice(S4, output= list)):
subsindets[3](Gp:-DrawSubgroupLattice(S4), specfunc(TEXT), k-> L[k], 2, applyop);

@jud I was using Maple 2019. I ran your worksheet in Maple 2019, and I get the correct graph. So I guess that there's a difference in Maple 2021, perhaps in the way that the numbers are represennted in the plot structure. I can't check Maple 2021 at the moment. However, if you execute the following commands and post the output, I may be able to diagnose the problem:

restart:
Gp:= GroupTheory:
S4:= Gp:-SymmetricGroup(4):
indets(Gp:-DrawSubgroupLattice(S4), specfunc(TEXT));

The output should look something like this:

{TEXT([-1.80739222823013, 1.], "2", FONT("HELVETICA", "DEFAULT", 8), COLOUR(RGB, 0., 0., 0.)), TEXT([-1.73719807243076, 2.], "15", FONT("HELVETICA", "DEFAULT", 8), COLOUR(RGB, 0., 0., 0.)), TEXT([-1.54919333848297, 1.], "3", FONT("HELVETICA", "DEFAULT", 8), COLOUR(RGB, 0., 0., 0.)), TEXT([-1.46993683051833, 2.], "16", FONT("HELVETICA", "DEFAULT", 8), COLOUR(RGB, 0., 0., 0.)), TEXT([-1.29099444873581, 1.], "4", FONT("HELVETICA", "DEFAULT", 8), COLOUR(RGB, 0., 0., 0.)), TEXT([-1.20267558860591, 2.], "17", FONT("HELVETICA", "DEFAULT", 8), COLOUR(RGB, 0., 0., 0.)), TEXT([-1.03279555898864, 1.], "5", FONT("HELVETICA", "DEFAULT", 8), COLOUR(RGB, 0., 0., 0.)), TEXT([-0.894427190999916, 3.], "26", FONT("HELVETICA", "DEFAULT", 8), COLOUR(RGB, 0., 0., 0.)), TEXT([-0.774596669241483, 1.], "6", FONT("HELVETICA", "DEFAULT", 8), COLOUR(RGB, 0., 0., 0.)), TEXT([-0.668153104781061, 2.], "18", FONT("HELVETICA", "DEFAULT", 8), COLOUR(RGB, 0., 0., 0.)), TEXT([-0.516397779494322, 1.], "7", FONT("HELVETICA", "DEFAULT", 8), COLOUR(RGB, 0., 0., 0.)), TEXT([-0.447213595499958, 3.], "27", FONT("HELVETICA", "DEFAULT", 8), COLOUR(RGB, 0., 0., 0.)), TEXT([-0.400891862868637, 2.], "19", FONT("HELVETICA", "DEFAULT", 8), COLOUR(RGB, 0., 0., 0.)), TEXT([-0.133630620956212, 2.], "20", FONT("HELVETICA", "DEFAULT", 8), COLOUR(RGB, 0., 0., 0.)), TEXT([0., 0.], "1", FONT("HELVETICA", "DEFAULT", 8), COLOUR(RGB, 0., 0., 0.)), TEXT([0., 1.], "8", FONT("HELVETICA", "DEFAULT", 8), COLOUR(RGB, 0., 0., 0.)), TEXT([0., 3.], "28", FONT("HELVETICA", "DEFAULT", 8), COLOUR(RGB, 0., 0., 0.)), TEXT([0., 4.], "30", FONT("HELVETICA", "DEFAULT", 8), COLOUR(RGB, 0., 0., 0.)), TEXT([0.258198889747161, 1.], "9", FONT("HELVETICA", "DEFAULT", 8), COLOUR(RGB, 0., 0., 0.)), TEXT([0.400891862868637, 2.], "21", FONT("HELVETICA", "DEFAULT", 8), COLOUR(RGB, 0., 0., 0.)), TEXT([0.516397779494322, 1.], "10", FONT("HELVETICA", "DEFAULT", 8), COLOUR(RGB, 0., 0., 0.)), TEXT([0.894427190999916, 3.], "29", FONT("HELVETICA", "DEFAULT", 8), COLOUR(RGB, 0., 0., 0.)), TEXT([0.935414346693485, 2.], "22", FONT("HELVETICA", "DEFAULT", 8), COLOUR(RGB, 0., 0., 0.)), TEXT([1.03279555898864, 1.], "11", FONT("HELVETICA", "DEFAULT", 8), COLOUR(RGB, 0., 0., 0.)), TEXT([1.20267558860591, 2.], "23", FONT("HELVETICA", "DEFAULT", 8), COLOUR(RGB, 0., 0., 0.)), TEXT([1.29099444873581, 1.], "12", FONT("HELVETICA", "DEFAULT", 8), COLOUR(RGB, 0., 0., 0.)), TEXT([1.46993683051833, 2.], "24", FONT("HELVETICA", "DEFAULT", 8), COLOUR(RGB, 0., 0., 0.)), TEXT([1.54919333848297, 1.], "13", FONT("HELVETICA", "DEFAULT", 8), COLOUR(RGB, 0., 0., 0.)), TEXT([1.73719807243076, 2.], "25", FONT("HELVETICA", "DEFAULT", 8), COLOUR(RGB, 0., 0., 0.)), TEXT([1.80739222823013, 1.], "14", FONT("HELVETICA", "DEFAULT", 8), COLOUR(RGB, 0., 0., 0.))}

@Ronan All of the commands that you've shown are top-level, not part of LinearAlgebra or any other package.

Please use the green up arrow on the toolbar to upload your worksheet. Copy and paste isn't really human readable when used directly on Maple output. 

And do we really need several pages of author acknowledgements at the beginning? 

@LegendCero 

Re question 4:

To get any solution, it's mathematically necessary that there be two more equations. They don't necessarily need to be differential equations. Any 2 equations (independent from the others) that use the same function variables, or some subset thereof, would be mathematically sufficient.

From looking at the physics and geometry of the situation, I think that the last two equations should look similar to the first two. Here's why: Consider the first equation: It describes the interactions between x_1 and its immediate neighbors---which are x_2 and y_1. Likewise, the second equation describes the interactions between y_1 and its immediate neighbors---x_1 and y_2. But the third and fourth equations describe the interactions between x_k or y_k (respectively) and their three (each) immediate neighbors---x_(k-1), x_(k+1), and y_k or y_(k-1), y_(k+1), and x_k (respectively). But at the right-side column, x_50 and y_50 only have two neighbors each. Thus, I think that the last two equations should be

m*diff(x_50(t), t$2) + k_g*x_50(t) = k_c*(x_49(t) - x_50(t))^3 + k_e*(y_50(t) - x_50(t)),
m*diff(y_50(t), t$2) + k_g*y_50(t) = k_c*(y_49(t) - y_50(t))^3 + k_e*(x_50(t) - y_50(t))

What do you think about that?

Re question 3:

If we assume that the springs have no transferable potential energy at time 0, the then total energy at any future time will be the kinetic energy of x_1 from its initial velocity. We'll divide all the energies by that to get a color to plot.

Re question 5:

I think that each mass's deflection should be measured relative to its own initial position. Thus, no step size h should be used, and 0 should be used for all initial positions.

If you agree with all of this, then I may now have enough information to get a numeric solution and plots.

I think that I can help you if you can supply the following information:

1. In terms of the variables (and parameters) that you used in the system of equations, what is the horizontal axis "Oscillator index"?
2. In terms of the variables (and parameters) that you used in the system of equations, what is on the color-scale legend bar that goes from 0.0 to 0.8? (I'd guess it's what you called "normalized energy".)
3. Does "normalized" mean that they're divided by the maximum value? And, if so, is it the maximum value within the boundaries of the plot? or some theoretical maximum? 
4. There are 2*n+2 function variables x_1, ..., x_(n+1), y_1, ..., y_(n+1), but only 2*n equations. You need two more equations.
5. You need 4*n+4 initial conditions such as n+1 initial positions (both x(0) and y(0)) and n+1 initial velocities (both x'(0) and y'(0)).
6. Does the v0 = 0.06 m/s mentioned in the plot titles refer to sqrt(x'(0)^2 + y'(0)^2)?
7. Is the "energy" that you're interested in kinetic energy, thus energy = m*v^2/2 = m*(x'(t)^2 + y'(t)^2)/2?
8. What specific value of n do you want?

@BekkiRR In the final plot command, find the line that begins contours= ..., and replace it with

contours= [$0..4] * 25000,

@Christian Wolinski From reading many of your previous Answers, I suspect that you are using a very old version of Maple.

You wrote:

• I would suspect the reason is you are not using diff and mtaylor does.

Of course we have found out in this thread that diff is much faster than D, but series is even faster. In your Maple, it may be true that mtaylor uses diff, but that isn't true in modern Maple; it uses series. Other than checking for several possible error conditions, handling the rarely used options for variable weights and returnorder, and doing some simplifications of the output, mtaylor is equivalent to this very short procedure:

BySeries:= proc(f::algebraic, V::list(name= algebraic), Ord::posint:= Order)
local t, v:= lhs~(V), v0:= rhs~(V);
    subs(O(1)= 0, t= 1, v=~ v-~v0, series(subs(v=~ v*~t +~ v0, f), t, Ord+1))
end proc
:

CodeTools:-Usage(BySeries(1/(2+x*y^2), [x=0, y=0], 30));
memory used=19.97KiB, alloc change=0 bytes, cpu time=0ns, real time=3.00ms, gc time=0ns

Is there a reason that you were suspicious about the correctness of the result? If so, I could likely explain it to you. Also, do you know how to get this solution by hand, starting with "separation of variables"?

@ecterrab I agree with this opinion that you just expressed:

• Add to that my opinion, that returning piecewise functions all around is not a convenient idea. Instead, I prefer presenting results with a convenient default and leave to the user to use assumptions for the other cases. 

Regarding those "other cases": Several high-level symbolic commands have an option for piecewise output for when you need all the detail. For example, solve has the parametric option (for polynomial systems with rational coefficients), and int has the allsolutions option (for definite integrals with parameters). I support the development of options like these. Although the amount of detail in these piecewise solutions is often overwhelming (especially in the solve case), a user can prune them to a manageable size programmatically or with assumptions.

This is a correction to @Christian Wolinski 's procedure:

P:= (x,y,x0,y0,N)->
local p,q,n,k;
    add(
        1/factorial(n)
        * add(
              binomial(n,k)*subs([p= x0, q= y0], diff(f(p,q),[p$(n-k),q$k]))
              * (x-x0)^(n-k)*(y-y0)^k, 
              k= 0..n
          ), 
    n= 0..N
):

@Anthrazit Is 11 larger or smaller than the problematic size?

I've occasionally had a similar problem with ordinary output (not in a container) which was corrected by increasing the zoom factor.

I see that the Mathematica output shows small plots of all the graphs. Something similar can be done in Maple like this:

PlotAll:= proc(L::list(Graph))
uses GT= GraphTheory, Pl= plots;
local n:= nops(L), c:= isqrt(n), r:= c;
    while r*c < n do r++ od;
    Pl:-display(
        Matrix(
            (r,c), 
            (i,j)-> 
                local k:= (i-1)*c+j; 
               `if`(k > n, plot(axes= none), GT:-DrawGraph(L[k]))
        )
    )
end proc
:
PlotAll(L2);

I don't know how to transfer this "array of plots" output to MaplePrimes, but you'll see it in a worksheet.