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.

I was afraid that might happen. Unfortunately, the syntax of the output option for Graph6 is not documented in the help pages, and I was making a guess. But this should work:

L2:= Import~(
    StringTools:-Split(ssystem("D:/nauty27r3/geng -c -b 6 -g")[2], "\n"), 
    source= direct, format= "Graph6"
); 
L2 := [`Graph 2: an undirected unweighted graph with 6 vertices and 5 edge(s)`,
 `Graph 3: an undirected unweighted graph with 6 vertices and 5 edge(s)`,
 `Graph 4: an undirected unweighted graph with 6 vertices and 5 edge(s)`, 
 `Graph 5: an undirected unweighted graph with 6 vertices and 6 edge(s)`, 
 `Graph 6: an undirected unweighted graph with 6 vertices and 5 edge(s)`, 
 `Graph 7: an undirected unweighted graph with 6 vertices and 6 edge(s)`,
 `Graph 8: an undirected unweighted graph with 6 vertices and 7 edge(s)`,
 `Graph 9: an undirected unweighted graph with 6 vertices and 8 edge(s)`,
 `Graph 10: an undirected unweighted graph with 6 vertices and 5 edge(s)`,
 `Graph 11: an undirected unweighted graph with 6 vertices and 6 edge(s)`,
 `Graph 12: an undirected unweighted graph with 6 vertices and 5 edge(s)`,
 `Graph 13: an undirected unweighted graph with 6 vertices and 6 edge(s)`,
 `Graph 14: an undirected unweighted graph with 6 vertices and 7 edge(s)`,
 `Graph 15: an undirected unweighted graph with 6 vertices and 6 edge(s)`,
 `Graph 16: an undirected unweighted graph with 6 vertices and 7 edge(s)`,
 `Graph 17: an undirected unweighted graph with 6 vertices and 8 edge(s)`,
 `Graph 18: an undirected unweighted graph with 6 vertices and 9 edge(s)`]

It definitely works in this case. It'll work in all cases (of format= "Graph6") if the newline character "\n" is always the separator between graphs and never occurs within an individual graph.

The strings that you're passing are just shell commands (for a Cygwin shell I assume). The Maple command for passing commands to a shell is ssystem. It will pass to the Windows/DOS Command/CMD shell. I doubt that you should begin the string with "!".

I agree with Tom Leslie that extending his printf loop to columns beyond the first differences would be difficult, that is, difficult if we're required to both generate and print the data in the same loop. But if you generate and store all the data before printing it, then it's a trivial one-liner. Add this to the procedures that I gave above:

PrintfTab:= (T::Matrix)-> printf("%9s\n", subsindets[2](T, numeric, "%2.6f", sprintf)): 

Example usage:

X:= <[$1..7]>:
PrintfTab(DisplayTab(DiffTab(1/~X), X));
 1.000000  1.000000                                                            
                    -0.500000                                                  
 2.000000  0.500000            0.333333                                        
                    -0.166667           -0.250000                              
 3.000000  0.333333            0.083333            0.200000                    
                    -0.083333           -0.050000           -0.166667          
 4.000000  0.250000            0.033333            0.033333            0.142857
                    -0.050000           -0.016667           -0.023810          
 5.000000  0.200000            0.016667            0.009524                    
                    -0.033333           -0.007143                              
 6.000000  0.166667            0.009524                                        
                    -0.023810                                                  
 7.000000  0.142857          

@abiodunowoyemi9 I don't know (mathematically) how to handle intial conditions for the derivative(s) of fractional DEs. Do you know?

Your labels say that the function has negative values in all four quadrants. That can't be true for the given function, and it's also contradicted by the colors of the regions.

Anyway, do you know the function values for which you want labeled contours?