Your object is not a Matrix but an Array.
Guessing this Array is N x N x P where P is, or not equal to N, you can use this
 

Download Array.mw

Look also to the ArrayTools package

Example

restart;
plots:-display(
  plot3d(x*exp(-x^2-y^2), x=-3..3, y=-2..2, color=red, style=surface),
  plot3d(exp(-x^2-y^2), x=-1..1, y=-1..2, color=blue, style=wireframe)
)

Download example.mw

As your procedure seems to be an iterative solver to find the value of TMP4 which, in some sense, is closed to the target TMP1, 
I suggest you to write the condition this way

 while abs(TMP4/(TMP1)-1) > 0.001 do

Then 

f(230, 34);

                          1.654689569

First point : in Maple the Hermite polynomials are mutually orthogonal with respect to the measure exp(-y^2) dy; (which differs from their definition in Statistics where there are mutually orthogonal with respect to the measure df(y), where f denotes the PDF of a standard gaussian random variable,  see for instance Hermite_polynomials).

Accordingly your definition of HN is wrong, and even doubly wrong because of the unjustified presence of an extra sqrt.
You should have define HN this way

HN := (n, y) -> exp(-y^2)*simplify(HermiteH(n, y)) / (sqrt(Pi)*2^n*n!):

Second point: in a reply @acer suggested you to use  [fsolve(simplify(HermiteH(7, x)))];  instead of    evalf(allvalues(RootOf(HermiteH(7, x), x))); and you answered him  "One solution is not enough, that's why I tried with solve. I need to get all the possible solutions, then discard complex ones and study the behavior of the ones that only provide pure real coefficients for the expansion".
It's pure nonsense!
All the zeros of Hermite polynomials are real. For a quick and simple proof look here the-roots-of-hermite-polynomials-are-all-real

What you want to compute are the values named "Gauss points", or "Gauss integration points", which are indeed the zeros of the Hermite polynomials..Have you ever heard of a Gauss quadrature method where the integration points would be complex ?
These zeros have symmetric values, implying that Hermite polynomials of odd degrees have 0 as a root.

You should follow acer's advice, all the more that the way you compute the Gauss points gives wrong results:

r := evalf(allvalues(RootOf(HermiteH(7, x), x)));
         0.8162878829, -1.673551629, -0.8162878829, 0.

Only 4 roots instead of 7!
(which explains why @tomleslie's solution contains so many undeterminate coefficients)

Here are the "corrrect" values

fsolve( simplify(HermiteH(7, x), 'HermiteH') );
 -2.651961357, -1.673551629, -0.8162878829, 0., 0.8162878829, 1.673551629, 2.651961357

Another way to get them is 

evalf([solve(simplify(HermiteH(7, x), 'HermiteH'))]);
sort(Re~(%));
            [                                -11    
            [0., 2.651961356 + 5.356554700 10    I, 

                                           -11    
              -2.651961356 - 5.356554700 10    I, 

                                           -10    
              0.8162878825 + 2.479096188 10    I, 

                                            -10    
              -0.8162878825 - 2.479096188 10    I, 

                                          -10    
              1.673551629 - 3.850404915 10    I, 

                                           -10  ]
              -1.673551629 + 3.850404915 10    I]
 [-2.651961356, -1.673551629, -0.8162878825, 0., 0.8162878825, 1.673551629, 2.651961356]

Finally: Here is a full computation of the ansatz u7
 

Download HermiteH_.mw

Oldies but Goodies: jresv48n2p111_A1b.pdf 

@perr7 

For m=0 the rhs should be exp(-x).
(So I repeat what I said in my initial answer: where does this equality come from?)

Download Equality_False.mw

Remark: in case you would think that Maple doesn't compute correctly LaguerreL polynomials, 
here is an alternative way to do the job.
It's based on the Rodrigues' formula for generalized Laguerre polynomials:

# Rodrigues' formula for Generalized Laguerre Polynomials 
# (aka LaguerreL(m, alpha, x))

GLP :=(m, alpha) -> piecewise(
                      m=0, 1, 
                      m::posint, exp(x)*x^(-alpha)/m! * diff(exp(-x)*x^(m+alpha), x$m), 
                      undefined
                    ):


# The two first GLP
eval(GLP(0, n-0));
simplify(eval(GLP(1, n-1)));
                               1
                             n - x
GLP_formula := m -> Sum(m!/n! * x^(m-n) * GLP(m, n-m)^2, n=0..+infinity):
for i from 0 to 2 do
  value(GLP_formula(i)):
  numer(%) / expand(denom(%)):
  print(i, simplify(%, hypergeom)): # simplification required for i >=3
end do:

# same results that those obtained using LaguerreL

A Generalized Laguerre Polynomial L(m, a, x)  (in LaTex notation $L_m^a(x)$ is defined only if a > -1.
In your formula a=n-m  which suggest that n should vary from m-1 (excluded ) to +oo for  L(m, a, x)  to exist (and not from 0 to +oo as you wrote):
Look to this:

restart:
with(orthopoly):
simplify(LaguerreL(2, 3, x), 'LaguerreL')
                                   1  2
                        10 - 5 x + - x 
                                   2   
simplify(LaguerreL(2, -0.999, x), 'LaguerreL')
                                                       2
       0.0005005000000 - 1.001000000 x + 0.5000000000 x 
simplify(LaguerreL(2, -1, x), 'LaguerreL')
                      LaguerreL(2, -1, x)

More of this the lhs of your formula can't be equal to exp(x) for it is a function of x and m.

Last but not least, the formula (in LaTeX style) seems wrong for x=0:

S := (m, N, x) -> add(m!/n!*x^(m-n)*LaguerreL(m, m-n, 0)^2, n=0..N):
p := 10:
S(p, p-1, 0);
S(p, p, 0);
                               0
                                       2
                       LaguerreL(10, 0) 

Could you be more explicit about your formula and give us some reference?

several possibilities

restart
with(Physics):
convert([a,b],`*`)
Error, unrecognized conversion

# force the use of the "default" arithop `*`:
convert([a,b],:-`*`)
                              a b

# or use a workaround
`*`(a,b):
mul([a,b]):

As the zeros of the ChebyshevT polynomial of degree N are 

zeros := sort([seq(cos((2*k-1)*Pi/2/N), k=1..N)])

it becomes interesting to find out how far it is possible to push the formal calculation.

The result is in the attached file.
The expression of the L1 norm is quite lengthy (uncomment the line K := ... to wtach it).


 

Download ChebyshevT.mw

I don't know what version of Maple you use but the integration can be done with Maple 2015.2 (look to the attached file).
Unfortunately alpha appears as a very complex function (MeijerG) of x_0 and there is no possibility in "invert" this expression to get x__0 as a function of alpha.

But this is not necesseray if you only want the graph of x__0 versus alpha.
Finally you could use CurveFitting:-Spline or Statistics:-LinearFit to buid an simplified representation of the "reverse" curve.

proposal.mw

About thirty years ago a colleague of mine did a PhD  titled "Formal tools for modeling in mechanics"
download (in french)
All the PhD is centered around the modeling and simulation, via Maple, of mechanical systems described by d'Alembert's, Lagrange's or Hamilton's equations.
Annex A contains a Maple V (five)  code dedicated to staellite motion.

We have lost contact since then, but tracking down his publications on the web could bring you some answers in english :-) )

Among other resources :
mapleprimes application center view.aspx

You will find two examples in the attached file: the first one sequentially displays random points on the unit sphere, the second one plots an heliw over the unit sphere.
Both of them are based on the use of plots:-animate and Statistics:-ScatterPlot3D (but this latter coud be replaced by pointplot for instance).
You said you have 5681 points and thus you have to ask yourself if you want to display them one at a time, or several at the time.
On the same way you will have to choose the number of frames you wnt to use.

In my examples I work with a smaller number of points (less than 500) and I decided to use a number of frames equal to the number of points. Each frame displays a new point.

Watch out: the mw file can be very big if you use a large number of frames (which is a reason why I loaded the link alone)

Points_on_a_sphere.mw

You can go a little bit further by doing this Lyapounov.mw  (iter=13, 500 pts in 84s).
It's still far from what is needed to obtain the blue curve in your post.



I guess the Iterator package (Maple >= 2019 (?),  I use Maple 2015 right now) could do the computations in a faster way.

Your error comes from the fact that you do not impose 4 ICs.
To plot y[3](tau) versus y[1](tau) add the option scene:

 

sys := [odey[1], odey[2], odey[3], odey[4]];

# Specify 4 ICs:

IC  := [ [y[1](0) = 0.1, y[2](0) = 0, y[3](0) = 0, y[4](0) = 0] ];

DEplot(sys, [y[1](tau), y[2](tau), y[3](tau), y[4](tau)], tau = 0 .. 2, IC, stepsize = .1, linecolor = blue, thickness = 2, arrows = medium, tickmarks = [0, 0], scene=[y[1](tau),y[3](tau)])

Or, fo several lists  of ICs
 

sys := [odey[1], odey[2], odey[3], odey[4]];
IC  := [ seq( [y[1](0) = 0.1*i, y[2](0) = 0, y[3](0) = 0.1*i, y[4](0) = 0], i=1..5) ];


DEplot(sys, [y[1](tau), y[2](tau), y[3](tau), y[4](tau)], tau = 0 .. 2, IC, stepsize = .1, linecolor = [red, gold, green, cyan, blue], thickness = 2, arrows = medium, tickmarks = [0, 0], scene=[y[1](tau),y[3](tau)])

Another way to do this, a little bit longer to use but which offers more possibility (IMO) , for instance in terms of adding a legend is this one (do not load DEtools):
 

sol := dsolve({sys[], y[1](0) = a, y[2](0) = b, y[3](0) = c, y[4](0) = d}, numeric, parameters=[a, b, c, d]):

and :

r := rand(0. .. 1.):

AllPlots := NULL:
for i in [seq(0..1, 0.5)] do
  for j in [seq(0..1, 0.5)] do
    sol(parameters=[i, j, i, j]):
    col := ColorTools:-Color([r(), r(), r()]):
    AllPlots := AllPlots, plots:-odeplot(sol, [y[1](tau), y[3](tau)], tau=0..2, color=col, legend=typeset([i,j]))
  end do:
end do:

plots:-display(AllPlots, gridlines=true, legendstyle=[location=right], size=[500, 400])

I used your worksheet and just add this at the end of it:

# you wrote "I can always define a function that computes the commutator"
# Yes, here is a naive example where neither the type of Phi nor e is checked

commut := proc(Phi, e)
  collect(simplify(mult(Phi, e)-mult(e, Phi)), [Dx])
end proc:

# your example
f := x -> sin(x)*cos(x)^8:
commut(A(s), f(x))

        /        2   \       7
        \9 cos(x) - 8/ cos(x) 



# "What command allows me to act with 'A(s)' on 'f' and get the resulting function?"
#
# I suggest you doing this

eval(diffop2de(A(s),u(x)), u(x)=f(x));
                           8         9           2       7
     s cot(x) sin(x) cos(x)  + cos(x)  - 8 sin(x)  cos(x) 

Example_mmdara.mw

@jennierubyjane  

Here is the answer to your removed question  (see @Preben Alsholm above)  
Comments and advices are in Maple 1D font within your worksheet named  ODEprobz-2 

Download ODEprobz-2_mmcdara.mw