The minimum shown is, of course, a bug. The correct value can be obtained by using the location option:

Min:= minimize(G(x,y),x=0..1,y=0..1, location);
                        9  
   Min := 5.838093910 10 , {[{x = 0.8949243201, y = 0.9715796788}, -0.7102080449]}

eval(G(x,y), Min[2][1][1]);
                          -0.710208045

Now, the erroneous value is still shown as the first number, but the true minimum is shown as the last number. Its correctness can be confirmed by direct evaluation and your plot.

This is the best that I could come up with so far. Its main drawback, to my mind, is that it creates a new DataFrame rather than just making a small change to the existing DataFrame. (In computer science lingo, we'd say that this doesn't operate inplace.)

data_2:= DataFrame(data_2, rows= [$1..upperbound(data_2, 1)])

My Answer is slightly different than @vv's because I assumed that you'd want to handle the more-general case of any number of rows in the new DataFrame. upperbound(..., 1) returns the number of rows, and the [$1.. ...] creates the list [1, 2, ..., nr], where nr is the number of rows.

I could easily have missed something easier. Although there are many help pages for DataFrame, they are not well organized IMO. In particular, the overview page ?DataFrame should at least list the name of every command specific to DataFrames, like the overview pages for most packages do.

Here's a way to find discrete local maxima. Suppose we have a list of (x,y)-points, sorted with respect to x, and we want the local maxima of the y-values.

#Generate test data by extracting a point-data matrix from a plot:
plot(sin(x), x= 0..5*Pi):
A:= indets(%, rtable)[]:

#Filter out consecutive equal y-values:
A1:= A[[seq](`if`(A[i,2]<>A[i+1,2], i, NULL), i= 1..upperbound(A)[1]-1)]:

#Find discrete local maxima:
convert(
    A1[[seq](
        `if`(A1[i,2] > max(A1[[i-1,i+1],2], i, NULL),
        i= 2..upperbound(A1)[1]-1
    )],
    listlist
);
        [[1.59043127800000, 1.], [7.87361657400000, 1.], [14.1568018700000, 1.]]

Download PyramidAngle.mw

Just by blind guessing, I added a sqrt(x) term to the argument of the exponential, and I got a near-perfect fit! Both this new model and your original can be linearized, which gives more-accurate results and which can be probabilistically quantified.

MaplePrimes won't let me display my worksheet, but you can download it.

Download BatteryDrainage.mw

Here's the plaintext code from the worksheet:

restart:
# Amps
A := <444, 218, 74, 312/10, 209/10, 119/10, 625/100>:
## Time until battery discharges to 10.5 volts`
T := <5, 15, 60, 3*60, 5*60, 10*60, 20*60>:
f := x-> a*exp(1-x/b+c*sqrt(x));

#The above model f(x) can be linearized:
ln(f(x)) = (ln(a)+1) + (-1/b)*x + c*sqrt(x);
LinearParms:= <ln(a)+1, -1/b, c>:
P:= Statistics:-LinearFit(
    [1, x, sqrt(x)], A, ln~(T), x, output= parametervector, summarize
);

Summary:
----------------
Model: 8.2359404+.12972487e-1*x-.58328210*x^(1/2)
----------------
Coefficients:
              Estimate  Std. Error  t-value   P(>|t|)
Parameter 1    8.2359    0.2723      30.2450   0.0000
Parameter 2    0.0130    0.0027      4.7533    0.0089
Parameter 3   -0.5833    0.0647     -9.0137   0.0008
----------------
R-squared: 0.9920, Adjusted R-squared: 0.9880

Model:= eval(f(x), solve({seq}(P =~ LinearParms), {a,b,c}));
plot([<A|T>, Model], x= 0..444); 

The piecewise-linear appearance of the plots can be removed by using odeplot with option refine. This will greatly reduce the number of points needed for smoothness, and it will figure out the number of variable-stepsize points on its own. It's an implementation of something akin to the adaptive plotting of the regular plot command.

All of your plots can be done and should be done with odeplot. The refine option must be used in conjunction with the range option in the dsolve command. It cannot be used with the numpoints option.

If you're going to do multiple plots over different ranges, then, of course, the range you pass to dsolve should be a superset of them all. I changed the order of your plots below to avoid an ugly-looking warning. The warning happens when a plot with a small range is followed by a plot with a larger range. The warning can be completely ignored if the larger range fits in the range passed to dsolve. I suspect that the warning message was put in the code before the refine option existed.

Several other tips:

• I relaxed the error tolerances to 1e-5 from your 1e-8. This is sufficient for accurate plotting, and it reduced the number of points in the largest plots by a factor of 4.
• Use the parameters option to dsolve. I've shown you this before.
• "Length of output exceeds ..." is not an error message or even a warning. It's an notification that Maple has sucessfully completed a calculation, but won't display it results because the amount of displayed data could seriously degrade the performance of the GUI display. The results of the calculation are still fully accessible.

MaplePrimes will not let me display this worksheet inline, I guess due to the large plots. You will see in the worksheet that they are all smooth.

Download OdeplotRefine.mw

Here is the complete code of the worksheet:

restart:
sys:= {
    diff(x(t),t) = x(t)*(r - b*x(t) - y(t)*(c + beta/(a+x(t)))),
    diff(y(t),t) = y(t)*(beta*x(t)/(a+x(t)) - mu)
}:
ics:= {(x,y)(0)=~ (0.2, 0.05)}:
Param:= [r,b,c,a,mu,beta]:
ParamV:= [1, 1, 0.01, 0.36, 0.4, 0.75]:
sol:= dsolve(
    sys union ics, {x,y}(t), parameters= Param, numeric, maxfun= 0,
    range= 0..20000, (abserr, relerr)=~ 1e-5
):
sol(parameters= ParamV);
plots:-odeplot(sol, [t, (x,y)(t)], t= 0..20000, refine= 1);
#Number of computed points in plot:
Npts:= P-> [rtable_dims]~([indets(P, rtable)[]]):
Npts(%%);
plots:-odeplot(
    sol, `[]`~(t, [x,y](t)), t= 0..400, axes= boxed, gridlines,
    color= ["Blue", "Red"], title= "Trajectories\n", legend= [x,y](t)
);
Npts(%);
eq_points_sym:= solve(eval(sys, diff= 0), {x,y}(t));
eq_points:= eval([eq_points_sym], Param=~ ParamV);
plots:-odeplot(
    sol, [x,y](t), t= 0..20000, refine= 1, color= red, thickness= 2,
    axes= boxed, gridlines, title= "Phase-space Diagram\n"
);
Npts(%);

#Task 2
ParamV:= [1, 1, 0.01, 0.36, 0.4, 1]:
sol(parameters= ParamV);
plots:-odeplot(sol, [t, (x,y)(t)], t= 0..20000, refine= 1);
Npts(%);
plots:-odeplot(
    sol, `[]`~(t, [x,y](t)), t= 0..400, axes= boxed, gridlines, refine= 1,
    color= ["Blue", "Red"], title= "Trajectories\n", legend= [x,y](t)
);
Npts(%);
plots:-odeplot(
    sol, [x,y](t), t= 0..20000, refine= 1, color= red, thickness= 0.9,
    axes= boxed, gridlines, title= "Phase-space Diagram\n"
);
Npts(%);
eq_points:= eval([eq_points_sym], Param=~ ParamV);

And here is my new version of the piecewisey plot that you showed in the Question. This is done with only 8733 points.

Here's how to get all the procedures neatly printed in the file:

SaveMyProcs:= proc(fname::string)
uses FT= FileTools, FTT= FileTools:-Text;
local 
    save_interf:= interface(
        longdelim= true, indentamount= 4, verboseproc= 2, screenwidth= 100
    ),
    fp
; 
    if FT:-IsOpen(fname) then FTT:-Close(fname) fi;
    try
        fp:= FTT:-Open(
            fname, create= true, overwrite= false, append= true
        );
        (P-> fprintf(fp, "%a:= %P;\n\n", P, eval(P)))~(
            select(type, [anames('alluser')], procedure)
        )
    catch: error
    finally 
        FTT:-Close(fp);
        interface(
            (longdelim, indentamount, verboseproc, screenwidth)=~
                save_interf
        )
    end try;
    return
end proc
:
SaveMyProcs("C://Users/carlj/desktop/MyProcs.mpl");

#The things in green are things that you might want to change.

#The text file looks like this:
MakeUnique:= proc(L::list, n::{nonnegint, identical(infinity)} := 1, f := x -> x)
    local e, T;
    option remember;
    T := table('sparse');
    [for e in L do if T[f(e, _rest)]++ < n then e; end if; end do];
end proc;

SaveMyProcs:= proc(fname::string)
    local save_interf, fp;
    save_interf :=
        interface(longdelim = true, indentamount = 4, verboseproc = 2, screenwidth = 100);
    if FileTools:-IsOpen(fname) then
        FileTools:-Text:-Close(fname);
    end if;
    try
        fp := FileTools:-Text:-Open(fname, create = true, overwrite = true, append = true);
        `~`[P -> fprintf(fp, "%a:= %P;\n\n", P, eval(P))](
            select(type, [anames('alluser')], procedure));
    catch:
        error ;
    finally
        FileTools:-Text:-Close(fp);
        interface((longdelim, indentamount, verboseproc, screenwidth) =~ save_interf);
    end try;
    return ;
end proc;

This yields 0:

simplify(residual) assuming (n, x, a, b, y(x)) >~ 0

Some positivity assumptions are usually needed to simplify expressions that have exponents with variables (n in this case---and note that n appears both as a base and an exponent here).

By the way, why do you use coulditbe for this? I don't think that you understand what it means: It's called existential quantification in mathematical logic. In other words, coulditbe returns true if there exists one or more specific numeric valuations of the variables which would make the statement true and satisfy the current assumptions. Isolated-point numeric evaluations are worthless for verifying a symbolic solution of a differential equation.

The command's name literally comes from the English phrase "Could it be true that..." which means exactly the same thing as "Is it possible that...", not "Must it be true that..." or "Is it necessarily true that...".

The source of your thread-safety problem and kernel crashes is not the ListTools procedures. ListTools:-Search is unequivocally threadsafe; ListTools:-MakeUnique is threadsafe in its single-argument form. There are many commands that are threadsafe that don't appear on the page ?index,threadsafe. It's quite sad that year after year passes, yet these don't get added to that list. For example, any expert reading the code for ListTools:-Search could tell you in 10 seconds or less that it is threadsafe. If you wish, you can replace these ListTools commands by

Search:= (x, L::list)-> local p; `if`(member(x,L,p), p, 0):
MakeUnique:= (L::list)-> (e-> ((thisproc(e):= ()), e))~(L):

I believe that the thread conflict in your makeproc comes from it returning its results via the global sparse Array j11 rather than by the usual return mechanism. 
 

As @dharr said, see lines 5-14 of `plots/odeplot`. There are numerous keywords that can be passed to a dsolve/numeric solution procedure instead of a number. One of these is "sysvars", and I think it shows what you want:

ds:= dsolve({diff(y(x),x)=y(x), y(0)=1}, numeric):
ds("sysvars");
                           [x, y(x)]

The issue that you're experiencing is specific to RootOf's use of _Z; it does not apply to _Z in general, and it has nothing to do with the underscore. Once you understand the issue, it's very easy to avoid.

First, let me explain bound variables and free variables. These are fundamental concepts of mathematical logic; they aren't specific to Maple or even to computer languages. For example, consider the commands 

int(a*x^2, x), int(a*x^2, x= 1..2),

in the absence of any surrounding context that may change what a and x mean. In both commands, a is a free variable and x is a bound variable. This is true regardless of whether they are local or global. The key thing is that x cannot be assigned any value (other than another name) without turning the command into nonsense.

Many commands require a bound variable or variables: int, diff, limit, sum, seq, solve (usually), plot (usually), RootOf, and many others. (In cases where the main argument is allowed to be a procedure, we can think of the procedure's parameters as the bound variables.) RootOf is unique among these (unfortunately) in the way that it handles its bound variable. If you do not specify a bound variable as the 2nd argument, it will look for global _Z in the first argument. If it's not there, it's an error. If you do specify a bound variable, RootOf will always change it to global _Z. Since _Z is protected, it should never be assigned a value, and you don't need to worry about that. So, in order for your example with local _Z to work, you need to put that _Z as the second argument of RootOf. (But there's really no reason to do this because you'd might as well just use global _Z, and, regardless, it's going to change whatever you use to global _Z anyway.)

The statement given at ?_ about those variables being reserved is ridiculously overstated. For one thing, it's not possible for a local variable beginning with _ to cause any problem.

You should never declare a global variable as global unless you intend to assign a value to it within the procedure containing the declaration.

You need to replace e^(-B.eta) with exp(-B*eta).

The number of boundary conditions required is the differential order, 6, plus the number of unspecified constants, e.

Maple does not "struggle with integer factorizations"; you just happened to stumble upon a particular number that exposes a bug. In the worksheet below, I randomly generate 1000 composites, each a product of a 9-digit prime and an 11-digit prime (like your example), and factor them all in 7 seconds.

Download ifactor.mw

When mod is used on a polynomial, such as t, it's mapped over the coefficients; so, t mod 10 becomes just t in the 2nd plot. This happens before plot even sees its arguments.

In the first plot, the special evaluation rule of seq prevents the simplification of t mod 10 to t.

This will do that job for any list:

ListTools:-Collect(L)