A solution is to ?overload the procedure `diff/Int`.

I coded the overload of `diff/Int` solution. It was actually much easier than I was anticipating, just 11 lines of code. This solution illustrates a very generally applicable technique, one that can be used to modify the behavior of almost any Maple library or even built-in procedure. This technique is not often discussed on MaplePrimes, so I suggest that most readers save this Answer in their Favourites.

restart;

`diff/Int/orig`:= eval(`diff/Int`):
unprotect(`diff/Int`);
`diff/Int`:= overload([
    proc(f::algebraic, x::{name, name=range}, M::(identical(method)=name))
        option overload;
        subsindets(
            diff(Int(f,x), args[-1]),
            specfunc(anything, 'Int'),
            subs([_DUMMY= args[3..-2]], J-> 'Int'(op(J), _DUMMY))
        )
    end proc,
    
    `diff/Int/orig`
]):

Test

I have only overloaded the case where the first argument to Int is type algebraic, the second argument is a name or name= range, and the third argument is method= name. All other cases should behave as they did originally. If you need the overload to work for more general cases of Int, please let me know.

Download diffInt_overload.mw

You wrote:

Notice that the order of the output from the EulerLagrange command when three functions are input does not correspond to the order of the functions listed, and when they are requested separately, the last function is not included:

The two phenomena are unrelated. As to the order: The EulerLagrange command always gives its result as a set, in curly braces. There is no way to control the order of elements in a Maple set: An internal lexicographic ordering is used so that set-membership checks can be done with binary search and set order is preserved across sessions. (Older Maple used an address-based ordering and set order varied by session.)

You wrote:

T3:=(1/2)*m*(diff(r(t),t)^2+r(t)^2*diff(theta(t),t)^2+diff(z(t),t)^2):
U3:=m*g*z(t):
L3:=T3-U3;

As for the second issue: Your equations above use diff(theta(t), t) but they never use theta(t) undifferentiated. This confuses EulerLagrange, and perhaps it should be rejected as invalid input. If you change diff(theta(t),t) to theta_prime(t) and do not include this in the solution variables, then the union of the individual results for z(t) and r(t) will be the same as the result when both of those are requested.

It can be done like this (I just picked a simple curve for the implicitplot).

P:= plots:-implicitplot(y= exp(x), x= -10..10, y= 0..10^4):
data:= op([1,1], P):
plots:-logplot(convert(data, listlist));
 

Change your final display command to

display([P, f(C)], viewpoint= ["circleright", frames= 32]);

See ?plot3d,viewpoint. The path of the "camera" can be arbitrary; I just picked an easy option above.

Add the option style= patchnogrid or style= wireframe to your implicitplot3d command.

It can be done, for example, like this:

# Returns an Array with n dimensions each with extent 1..m
DynArray:= (n,m)-> Array((1..m)$n, args[3..-1]):
A:= DynArray(3,4);
              

Be careful with this! A has m^n elements. 2^33 = 8 Gig of elements.

plot3d(unapply(S,x,y), 0..1, 0..1, labels= [x,y,``]);

To look at the contents of a large matrix, right-click on it and select Browse.

As to exporting to Excel: Are you saying that some entries in your Matrix are themselves lists? Do the lists have more than one element? If yes, do they all have the same number of elements?

This may be too tedious to be worth it: Use package ?Worksheet. This will read in a worksheet and convert it to a very large function tree. For the first case, find, by whatever ad hoc method, the function name that corresponds to the headings that you want. Then find, by experimenting with op, which operand or suboperand of that function is the heading string that you want.  Now you can extract all occurences that function name from any worksheet using indets(..., specfunc(anything, ...)), and then extract the heading string from those with op.

Here are two versions of contour plots of the lattice polynomial. To deal with the extremely large values, I took the 25th root of the polynomial (since the polynomial is degree 25). One problem with the Mathematica plot is the both extremely positive and extremely negative values are shown in white.

Contourplots.mw

P:= mul(mul((x+I*y)-(a+b*I), a= -2..2), b= -2..2):
plots:-contourplot(
   surd(Re(P),25), x= -3..3, y= -3..3,
   grid= [100,100], contours= 20,
   filled= true, coloring= [brown, green],
   axes= boxed
);

plot3d(
   surd(Re(P), 25), x= -3..3, y= -3..3,
   grid= [100,100],
   orientation= [-90,0,0],
   style= "patchnogrid", shading= zhue
);

LinearAlgebra:-IntersectionBasis is a great one-step solution to this problem. But for the sake of pedagogy, which may satisfy the Asker, and to reinforce my own learning, I present here several more solutions of this problem expressed in terms of more elementary operations.

Download IntersectionBasis.mw

Here is the smoothest plot that I can get without using image-processing tools (expanding on an idea of Markiyan's presented in the thread he refers to above).

P3:= plots:-implicitplot(
   y= sin(x), x= 0..2*Pi, y= -1.05..1.05,
   grid= [32,32], gridrefine= 4, crossingrefine= 4, resolution= 2^12,
   rangeasview, thickness= 2, color= brown,
   axis[1]= [tickmarks= piticks]
):
op([1,1,2,1,2], P3); #Check number of points used.
                               84
print(P3);

Assign the solution given by dsolve to a variable:

   soln:= dsolve(..., x(t));

where you have to substitute your actual dsolve command. Then, to plot it, do

   plot(eval(x(t), soln), t= 0..10);

If this doesn't work for you, please post or upload your code!

eval(V[dis], Units:-Unit= 1);

(A comment irrelevant to using the above command but essential to understanding why it works: The global Unit is initially assigned (i.e. before any packages are are loaded) to

    proc() Units:-Unit(args) end proc,

and `print/Unit` is initially assigned to

    proc() 
       local Unit;
       if IsWorksheetInterface('Standard') or _EnvInRecursiveUnitsModulePrint = true  then Unit(args)
       else Unit([args])
       end if
    end proc.

)

You need to make a copy of the names x, y, and r with convert(..., `local`). You store the copy under a different name, and when that name is invoked, it produces the copy of the original name, which has no value assigned, even if the original original was assigned. Whew! That's probably very confusing, so here's an example:

restart;
B:= convert(a, `local`);
                               a
a:= 1;
                               1
B;
                               a
a;
                               1

The second thing that's needed for your plot is that we want the point to display as (x₀, y₀) rather than as x₀y₀. So we invoke that with the empty-name function as ``(x[0], y[0]), and the quote marks do not appear in the typeset display. One other change that I made to your code is to put the alignments inside the brackets [] of the individual text elements, because I think that the point should be aligned {below, left} and that the radius should be aligned {above, right}.

Download plotcircle.mw

Note that the gridlines and subticks do not occur in the worksheet; that is just an effect of uploading to MaplePrimes.