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

Note that Maple's plotting commands produce somewhat high-level and human-readable data structures before they're displayed as plots. The relevent data can be extracted from those, as in the following example. The particular data that you asked about is in two-column matrices, which there may be several of.

P:= plot(x^2, x= -2..2); #Store plot in variable P.
op(P); #Show the entire data structure.
data:= indets(P, And(rtable, (posint, 2) &under upperbound)); #Extract the 2-column matrices

In this case, there's one such matrix, which I now extract from the set produced by the previous command:

Curve:= data[1]; #1st (and only) item in the set

Curve contains about 200 points (version dependent) in this case. Mine has exactly 200. Each row is an (x, y) pair in Cartesian coordinates (even if you used a different coordinate system, for example, polar).

Here's how to index matrices. The 1st entry in the 99th row is Curve[99, 1]; the 19th row is Curve[19]; the 2nd column (the y-coordinates) is Curve[.., 2]; etc. The are many indexing variations possible.

The command usually recommended for extracting data from a plot is plottools:-getdata. However, I strongly prefer the method that I described above.

If you make a numeric dsolve call such as

P:= dsolve({diff(y(x), x) = y(x), y(0)=1}, numeric, method= rkf45)

then P will be a procedure. List its contents with

showstat(P)

For a procedure of this size, showstat will give you more-human-readable results than print(P).

The procedure P has a **local** Array _dat, which has another procedure stored in _dat[1]. This latter procedure does the bulk of the work.

You said that this was in a procedure. In the procedure's declarations area, include

uses plots;

Then within that procedure, any reference to a command from package plots will not need the plots:- prefix.

I'd recommend against using alias or macro inside a procedure, but I'm only guessing that they might cause problems, potentially by altering the global state. I don't know from experience.

Regarding your statement odeplot:= plots:-odeplot: If odeplot is a local variable of the procedure, there would be no problem with doing it that way. But I think that uses is very, very slightly more efficient because there's no need to allocate and dereference a local.

@Ronan You asked:

• Q1. if A :=[2,3] how do I get the procedure to detect and print that it is _T2L?

See my additions to the attached worksheet/module, in green again. I also modified ModuleUnload so that you won't get those distracting but totally harmless warning messages about overwriting the types.

• Q2. why does it need [ModuleLoad() right before end module]?

Good question. It's just for editing, testiing, rewriting the code. A ModuleLoad is automatically called when the module is loaded from a library/archive. When we're rewriting/editing the code in a worksheet session, the new code is not in a library, so that automatic execution of ModuleLoad doesn't happen. But we need the effects of ModuleLoad() to test the code.

Download Module_Test_Types_2.mw

I guess that this is just a a very simplified example of what you actually want to do?

Matrix indices start at 1, not 0. So perhaps you want an Array or a table instead of a Matrix?

A Matrix can be indexed A[i, j] or A(i, j), but these mean slightly different things. A[i, j] assumes that i and j are valid indices. A(i,j):= ... will create them as valid indices if need be.

It's much easier to let the overload command do all the type-checking work rather than trying to do it in ModuleApply. Here is your module with my corrections. What I added is in green; what I removed is in orange.
 

Download Module_Test_Types.mw

The algorithm shown by @dharr is fairly well known; I remember reading it in a magazine (perhaps Games) as a child. The following property of that algorithm is much less well known. It doesn't help much when applying the algorithm by hand, but it allows for a big simplification in a computer implementation. It's this: The steps where one needs to go down instead of up and right due to a previously used position or going off the upper right corner happen iff the previously used value is a multiple of the order.

Also, for the verification procedure, I added a check that the elements are all different, which is usually considered part of the definition of a true magic square.
 

Download MagicSquare.mw

Like this:

ex:=
    G1*P3 + G1*P5 + G1*P6 + G2*P3 + G2*P6 + G3*P2 + G3*P5 + G4*P2 + G4*P3 + G4*P5 
    + G4*P6 + G5*P2 + G5*P3
    + G5*P5 + G5*P6 + G6*P3 + G6*P6 + G7*P2 + G7*P5 + G8*P2 + G8*P3 + G8*P5 + G8*P6
:
V:= indets(ex)[]:
codegen[optimize](unapply(ex, [V]), tryhard)(V);
    (G3 + G4 + G5 + G7 + G8)*P2 + (G1 + G3 + G4 + G5 + G7 + G8)*P5 + (P3 + P6)*(G1 + G2 + G4 + G5 + G6 + G8)

By the way, collect will also accept a list of variables: 

collect(ex, [P3,P6,P2,P5]);

It can be done with a one-line procedure:

Rose:= (n::integer)-> plots:-polarplot(sin(n*t), filled, color= `if`(n::odd, pink, aqua));

Test:

Rose(2);  Rose(3);

Given a planar graph G, the command

PD:= GraphTheory:-PlaneDual(G)

returns the graph PD whose vertices are the regions (a.k.a. "faces") of G and such that {u,v} is an edge of PD iff the corresponding faces of G share an edge. Unfortunately, this command doesn't give the face-to-vertex correspondence; however, it may be enough for your purposes.

The list of lists FL of vertices of G that define its faces can be obtained by 

GraphTheory:-IsPlanar(G, 'FL')

Vertices of degree 0, 1, or 2 in G can be a nuisance. I know this from long experience; I don't know whether they have any bad effect on these particular commands. Those vertices usually add nothing of interest mathematically to the study of planar graphs (or graphs embeddable on other surfaces), but they make many algorithms more complicated. I recommend that they be removed from G. Those of degree 0 or 1 can be removed outright. For a degree-2 vertex v, remove it and connect the two other vertices of its neighborhood if they aren't already connected. Continue recursively removing vertices of degree < 3 until there are no more.