Looking at the code with showstat(GraphTheory:-IsHamiltonian) clearly shows that there's no mention of tsp nor use of any algorithm akin to solving a traveling-salesman problem. So, the help page is clearly wrong. There is a procedure GraphTheory:-TravelingSalesman, it has a help page, and I tested it 6 days ago (in response to this Question: "Traveling-salesman problem"). 

It seems necessarily inefficient to me to use a minimizing algorithm if you just want any Hamiltonian circuit. However, if you want, I could easily (30 - 60 minutes) write a patch for you to implement method = tsp in IsHamiltonian. So, let me know if that's what you want.

The StringTools package can process so-called regular expressions (such as are commonly used in Linux shell scripts). Let's say that you want to extract a number (possibly including periods) in parentheses from the end of a string.

StringTools:-RegMatch("\\([0-9.]*\\)$", L1, 'c');
                              true

c:= c[2..-2]; #Strings can be indexed like Arrays and lists.
                          c := "2.13"
#Even more slick:
StringTools:-RegMatch("\\(([0-9.]*)\\)$", L1, 'c', 'c'), c;
                          true, "2.13"
​​​

All that you need is 

save u11, "C:/double/u11.txt";

The formatting, printing, and file closing statements are all unnecessary. You can also save all the u variables to one file like this

save u1||(1..10), "C:/double/u.txt";

All necessary punctuation (u11:= ..., semicolons, line breaks) will be included so that the file can both be read by a human as text and read back into Maple as executable code.

I've never a practical (useful) case of a nonsingular matrix with symbolic entries whose Moore-Penrose pseudoinverse could be calculated in a reasonable amount of time. 

The analog of "gradient" for a vector field is called jacobian. It is a matrix. It is available as command Jacobian in the VectorCalculus and Student:-VectorCalculus packages.

If the animation is produced with the Explore command, as shown below, the frames are displayed "on the fly" (i.e., as they are produced). So, unlike Maple's other animation commands, this has the benefit of starting play immediately. The memory used by the kernel for this is trivial, under 50 Mb; however, the memory used by the GUI grows unboundedly. I killed the below (just by pressing the stop button under the animation) when my GUI was using about 8 Gb. Still, it displayed a few thousand frames. (If you delete the Explore window after stopping it, those Gigs of memory will soon be garbage collected and returned to the O/S.)

The trick that I show in the code below---using ArrayTools:-Alias to create a direct memory link to the animation frame---can only work if the frames are displayed on the fly.

restart:
#
#Build unit sphere as static background for the animation. The back half is
#opaque, and the front half is wire mesh, so we can see inside.
#
Sphere:= plots:-display(
    plot3d~(
        1, [-Pi..0, 0..Pi], 0..Pi, coords= spherical, style=~ [patch, wireframe],
        lightmodel= light1, grid= [100$2], axes= normal, orientation= [60,80,10]
    )
):
N:= 5000: #number of points and animation frames
#Construct matrix each column of which is a random point on unit sphere:
Pts:= rtable(1..3, 1..N, frandom(-1..1), datatype= hfloat, subtype= Matrix):
Pts[3]:= csgn~(Pts[3])*~(1 -~ Pts[1]^~2 -~ Pts[2]^~2)^~(1/2):
Z:= Vector(3, 0, datatype= hfloat): #the origin
Frame:= plots:-display( #Construct one animation frame 
    Sphere, plots:-pointplot3d(<Z | Pts[.., 1]>, style= pointline, color= red)
):
#Create a direct memory link to the 2x3 matrix inside the pointplot3d structure.
#Thus, we can change the frame simply by changing the matrix.
A:= ArrayTools:-Alias(indets(Frame, Matrix)[]):
#
#Procedure that changes the frame:
P:= k-> 
    if k::posint then :-A[2]:= :-Pts[..,k]^+; :-Frame
    else 'procname'(k)
    fi
:
Explore(
    P(k), parameters= [[k= 1..N, shown= false, animate]], autorun, numframes= N,
    initialvalues= [k= 1]
);

One frame:

I suggest that you create an object (see help page ?object) to represent the 32-bit floats. The internal storage of the object should use strictly nonnegative integers. The floats can be packed into a 1-word integer. The Bits package can be used for reasonably efficient sub-word-sized manipulations. All arithmetic operations can overloaded to work with the floats.

How much programming experience do you have, and does that include any experience in object-oriented programming (OOP)?

The purpose of my suggestion is to facilitate experimentation as you develop the firmware code. Execution speed is not a concern.

If your document begins with restart (which is almost always a good idea), then the code for a procedure must come before the attempt to use the procedure.

The "unable to parse" error is caused by the semicolon at the end of the first line of myProc. Consider these three variations:

1. myProc:= proc(pA)
2. myProc:= proc(pA)::string;
3. myProc:= proc(pA);

The first two of those are valid syntax; the third is not (using 2D Input). 

Try rtable_eval(eval(solution, {t= 1})).

Not equals: <>

Solving an equation is almost equivalent to solving the numerator of the difference of the two sides. This gives a usable solution:

solve(numer(p0 - p1),  y);

It still gives that warning, but a warning is not an "error", and "may have been lost" is not necessarily "were lost". 

You should collect the data inside the loop into matrices, then use a plot command after the loop.

This does side-by-side plots for your 3 values of k1. In your data-collection loop, the lines that I added or modified are followed by ####. For brevity here, I omitted lines that weren't needed to get the plots; however, the code executes fine if those lines are included.

#Data-collection loop:
beta:= {$1..8}:  ####
for k1 in K do
    plotdata[k1]:= Matrix(nops(beta), 4); ####
    for i to nops(beta) do ####
        b:= beta[i]; ####
        Etemp := eval(
            `πcentral`(p, e, z), 
            [A = 0, B = 2, alpha = 50, beta = b, c = 5, k = k1, mu = 10, v = 1]
        );
        Soln:= NLPSolve(Etemp, p = 5 .. 500, e = 0.5 .. 10, z = 0 .. 10, maximize);
        (etemp, ptemp, ztemp):= eval([e, p, z], Soln[2])[]; ####
        `πtemp` := Soln[1];
        y_temp:= eval(y(p, e), [alpha = 50, beta = b, p = ptemp, k = k1, e = etemp]);
        q_temp:= y_temp + ztemp;
        plotdata[k1][i] := <b | ptemp | etemp | q_temp> ####
    od
od:
plots:-display(
    <seq(
        plot(
            [seq](plotdata[k1][..,[1,j]], j= 2..4), legend= ['p', 'e', 'q'],
            style= pointline, symbol= [box, diamond, solidcircle], 
            labels= ['beta', ``], title= sprintf("k1 = %a", k1), axes= boxed
        ),
        k1= K
    )>^%T
);

Yes, this is called a traveling-salesman problem. In the code below, the variable met indicates which norm (or metric) is used to measure the distance between points. You may not know that your points are already sorted with respect to the standard Euclidean (or "2") norm. If the 1-norm (sum of distances in each dimension) is used, the order is different. The code below takes about 4 minutes. This problem is a canonical example of an NP-complete problem, for which there's no known fast algorithm for guaranteed shortest path.

In my initial matrix below, note that I trimmed off the repetition of the 1st point as the last point.

met:= 1: #Use 1-norm.
LT:= Matrix([
    [1, 1, 0], [3, 2, 0], [2, 3, 1], [1, 3, 1], [1, 4, 2], [1, 5, 0], [2, 5, 1], 
    [3, 5, 3], [4, 5, 3], [3, 4, 5], [4, 3, 2], [5, 3, 2], [5, 4, 1], [5, 2, 1],
    [5, 1, 0], [4, 1, 2], [3, 1, 2], [2, 1, 2], [1, 2, 1]
]):
W:= Matrix( #"Weight" matrix, or pairwise distances
    upperbound(LT)[1]$2, (i,j)-> LinearAlgebra:-Norm(LT[i]-LT[j], met),
    shape= symmetric, datatype= hfloat
):
(d,P):= CodeTools:-Usage(GraphTheory:-TravelingSalesman(GraphTheory:-Graph(W)));
memory used=21.20GiB, alloc change=0 bytes, 
cpu time=3.87m, real time=3.60m, gc time=33.45s

d, P := 42., [1, 19, 4, 3, 7, 6, 5, 10, 8, 9, 11, 12, 13, 14, 15, 
  2, 16, 17, 18, 1]

#d = 42 is the total distance.

LT__sorted:= convert(<LT, LT[1]>[P], listlist);
LT__sorted := [
    [1, 1, 0], [1, 2, 1], [1, 3, 1], [2, 3, 1], [2, 5, 1], 
    [1, 5, 0], [1, 4, 2], [3, 4, 5], [3, 5, 3], [4, 5, 3], 
    [4, 3, 2], [5, 3, 2], [5, 4, 1], [5, 2, 1], [5, 1, 0], 
    [3, 2, 0], [4, 1, 2], [3, 1, 2], [2, 1, 2], [1, 1, 0]
]

There are cases where a discrete, plot-based technique such as you suggest is necessary, but this isn't one of them: Standard algebraic, analytic, numeric techniques will work for this. I think that you should use them when they work, because they'll give you much more mathematical flexibility. 

fsolve(1 - binomial(180000,r)*r!/180000^r = 0.25, r= 1..infinity);
                          322.2204929

This works because the right side of the equation is a continuous, increasing function of r for real r >= 1; r isn't restricted to integers.

There may be some cases where you'd need to use a plot-based technique (such as you suggest) to answer that question, but this isn't one of them. Here, straightforward algebraic, analytic, numerical techniques work, and I think it's more mafhematically robust for you to use them. 

fsolve(1 - binomial(180000,r)*r!/180000^r = 0.25, r= 1..infinity);
             322.2204929