Here is how you can get a package to automatically load all of its subpackages when it is loaded. My method is simpler than VV's in the long run because nothing needs to be changed as new subpackages are added:

restart;
OM:= module()
option package;
export
    IM:= module()
    option package;
    export add_three:= x-> x+3;
    end module,

    init:= proc($)
    local p;
        for p in exports(thismodule) do
            if thismodule[p]::':-package' then
                parse(sprintf("with(OM:-%a);", p), ':-statement')
            fi
        od;
        return
    end proc
; 
end module
:
with(OM);
                           [IM, init]
add_three(5);
                               8

Notes on that code:

1. The procedure must be named init, and it must be an export of the main package. If such a procedure exists, it'll be called automatically when package is loaded with with.

2. An init procedure is similar to, but not identical to, a ModuleLoad procedure. The former is called automatically when a package is loaded with with, while the latter is called automatically when a module (regardless of whether it's also a package) is read from a library.

3. The encapsulation of the with commands as strings executed by parse bypasses the usual prohibition against using with inside procedures and modules. Executing with parse is the same as executing at the top level.

4, My init requires no foreknowledge of the names of the subpackages nor of their exports. Indeed, it would be safe to keep it even if there were no subpackages.

5. My init can be used nearly verbatim in any package. The only thing that would need to be changed is the single occurence of the name of the outer package.

Here are two ways: nested piecewise or boolean combinations of the conditions:

f1:= piecewise(
    x<0, 0, x<1, piecewise(y<0, 0, y<1, x+y),
    x<2, piecewise(y<0, 0, y<1, x-y)
);

f2:= piecewise(
    Or(x<0, y<0), 0,
    And(x<1, y<1), x+y,
    And(x<2, y<1), x-y
);

Edit: Corrected direction of inequalities.

Here's an easy correction to your procedure:

LT:= (e::algebraic, t::algebraic)-> select(has, e + _x, t);

And, even better, make _x local. In Maple 2019 or later, that can be done in one line as

LT:= (e::algebraic, t::algebraic)-> local x; select(has, e+x, t);

You should keep in mind that has does a "deep" search. For example, 
has(y*sin(1/(1+f(t))), t) returns true. If you'd prefer a breadthwise-only search through the operands of the term, let me know. 

Yes, modules and procedures are first-class objects (see Wikipedia article) in Maple, meaning they can be placed anywhere any other expression can be placed. In contrast, libraries are not first-class objects.

Your usage of unapply is quite awkward, although I understand (from your Reply to mmcdara) why you did it. However, it's not needed: Create u as an expression, not a function, then subs for the X__4 or X__2 when that's syntactically required for integration or summation.

The symbolic integration and summation of all 800 terms is trivial, and can be done in 0.3 seconds:

u:= exp(-(X__4-Lmu)^2)*exp(-(X__2-Wmu)^2):
U:= CodeTools:-Usage(
        sum(int(subs([X__4= 5*k, X__2= x], u), x= -1..1), k= 1..800)
):
memory used=18.21MiB, alloc change=67.01MiB, 
cpu time=266.00ms, real time=264.00ms, gc time=15.62ms

The limits of integration could be made symbolic also. The entire integrated and summed expression above simplifies to 

sqrt(Pi)/2*(erf(Wmu+1)-erf(Wmu-1))*Sum(exp(-(Lmu-5*k)^2), k= 1..800);

Conversions between lists, sets, sequences, tables, and rtables are a constant requirement of Maple syntax; there's nothing special about the plotting commands that you mentioned. There's some such conversion in nearly every line of code that I post on MaplePrimes. Perhaps you'd feel better about it if you didn't use the lengthy syntax that you do. The conversion of an rtable R to a list requires only 7 additional characters:  [seq](R).

The code below creates and solves the same matrix-vector sys as you had. Together, that takes 1 line of code. Then it produces all 3 plots that you specified using [seq](sys) and the exact arguments that you used plus a title for each. The entire thing, including tab indentation and newlines, is 267 characters. (To contrast that, this paragraph has 362.)

restart:

V:= [x,y]; F:= <V(t)>; sys:= diff(F,t) - <1,2;0,3>.F; dsolve(sys);

for P in [dfieldplot, phaseportrait, DEplot] do
    (print @ DEtools[P])(
        [seq](sys), V(t), t= 0..4, 
        `if`(P=dfieldplot, [][], [V(0)=~ [1,0]]), V[]=~ -4..4, 
        title= sprintf("DEtools:-%a", P)
    )
od;

The command implicitdiff handles this. See its help page, ?implicitdiff.

A procedure for it:

MROUND:= proc(X::realcons, m::posint)
local x,r;
    if X<0 then return -thisproc(-X, m) fi;
    x:= trunc(X);
    r:= irem(x,m);
    if r + frac(X) < m/2 then x-r else x-r+m fi
end proc
:

This bug seems to be the result of a bit of arrogance on the part of the author of SolveTools:-CancelInverses. Look at lines 65-66. That person believes that there is an intermittent bug in simplify such that sometimes it returns expressions with a free _Z (which has been extracted from a RootOf). And, there may in fact be such a bug, but, if there is, it plays no part in the current example. The arrogant part is that that programmer does not check whether the expression had a free _Z before it was passed to simplify. Indeed, in this case, the expression has a free _Z even at the point that it's passed to CancelInverses. This makes it particularly easy to write a workaround for this bug:

restart:

interface(warnlevel=4):
kernelopts('assertlevel'=2):

CI__orig:= eval(SolveTools:-CancelInverses):
unprotect(SolveTools:-CancelInverses):
SolveTools:-CancelInverses:= overload([
    proc(e::satisfies(e-> depends(e, _Z)))
    option overload;
    local _A;
        eval(CI__orig(eval(e, _Z= _A), _rest), _A= _Z)
    end proc,

    CI__orig
]):
protect(SolveTools:-CancelInverses, CI__orig):

expr:= exp(x)*sin(y)-3*x^2+(exp(x)*cos(y)+1/3/y^(2/3))*Z = 0;
solve(expr, y);

The representation z <= sqrt(-x^2-y^2+9) of a hemisphere is in cartesian, not spherical, coordinates. So, if you remove coords= spherical, you'll get something better. But I think that this is even better:

plots:-display(
    plot3d(
        [3, theta, phi], theta= -Pi..Pi, phi= 0..Pi/2,
        coords= spherical
    ),
    plot3d(
        [r, theta, 0], r= 0..3, theta= -Pi..Pi, coords= cylindrical
    ),
    scaling= constrained
);

I think that's it's only really possible to plot surfaces, not solids. So, in the above, I "capped" the bottom of the hemisphere.

Yes, your data can be put into a matrix or array (if it isn't already so), and, like any Maple variable, that array can be saved to a file (with the save command) and retrieved in any other session (i.e., after a restart or even in a completely different worksheet) with the read command. If I do this with an array of the size that you describe, and I use some efficiency options (described below), then the save takes 1.3 seconds and the read takes 1.7 seconds. (I'm using an SSD drive, and those times are CPU-bound. I think that it'll take longer with a mechanical drive.) Here's an illustration:

restart:
#random array simulating your conditions:
A:= rtable(1..434000, 1..28, frandom(0..1), datatype= hfloat):
A[42, 23]; #some random entry for verification
                       0.0293574166199714

#Set a default directory for file I/O. You need to change this to
#something relevant to your computer:
currentdir("/Users/carlj/desktop"):

st:= time[real]():
save A, "mydata.m":
time[real]()-st;
                             1.234

restart:
A[42, 23]; #verify memory is cleared:
                           A[42, 23]

currentdir("/Users/carlj/desktop"):
st:= time[real]():
read "mydata.m":
time[real]()-st;
A[42, 23];
                             1.606
                       0.0293574166199714

The three verification steps are only there for educational purposes, and the four timing steps are only there for experimental purposes. None of those seven steps are actually needed.

Having this array in memory will not have any effect on the efficiency of Maple's GUI display. The array (like other ordinary variables) is not saved with the worksheet.

The two efficiency options mentioned above are

1. The file's name ends with ".m". This saves the data in Maple's internal format, essentially a character-encoded binary.
2. The array has a hardware datatype. I used hfloat, equivalent to float[8] or what's commonly called "double precision". Because of this, the array cannot contain the string data. That data can be saved in another array, and the two arrays can be stored in the same file. Or, the timestamp data could be encoded numerically (as its difference from a specific "base" time) and stored in the numeric array. The numeric encoding of timestamps is fairly easy in Maple.

Without these two efficiencies, the times are much longer (a little over a minute), which is still less than reading the CSV.

I'll assume that you realize that both results are correct, so I won't discuss that further (unless you ask).

First, the difference has nothing to do with using Interactive rather than Minimize. If the first case is entered as 

Minimize(x + y, {1 <= x + y}, x= 0..1, y= 0..1);

then you get the same result as your Interactive call.

Second, additional information about what's going on in the background can be obtained by setting

infolevel[Optimization]:= 5:

before making the call to Interactive, Minimize, or any other Optimization command. The extra information appears in the worksheet, not in the Maplet window (the dialog box that Interactive uses). One thing revealed by this information is that Maple knows that the solution is not unique; it refers to it as a "weak solution".

Setting a particular infolevel only needs to be done once per session (i.e., per restart).

Third, another interesting case which produces yet another solution is

Minimize(x + y, {x <= 1, y <= 1, 1 <= x+y}, assume= nonnegative);

The reason for these differences is that a simple inequality such as x <= 1 is treated just like a general multivariate linear inequality such as 1 <= x + y. Each such inequality leads to an extra row in the matrix. But there are more-efficient means of handling constraints entered as bounded ranges, such as x= 0..1; those don't lead to extra rows in the matrix. I recommend using this bounded-range form whenever possible, even for nonlinear and integer-linear problems.

The $ is the sequencing operator (see help page ?$). It still works in current Maple. The simple usage shown---a$n, where n is a nonnegative integer---produces a sequence of n copies of a. Thus [a$n] is a list of n copies of a. 

I'll assume that you're also removing usages of the old linalg package. You certainly should if you're not already. In that case, replace rowdim(A) with upperbound(A)[1] or LinearAlgebra:-RowDimension(A). Either will work; the first is more efficient.

The command plot3d is only for plotting surfaces. To plot curves in 3D, use plots:-spacecurve. I also added a slight transparency to the surface to make the ellipse fully visible.

f:= (x, y)-> x^2 + y^2 - 12*x + 16*y:
plots:-display(
    plot3d(f(x,y), x= -9..9, y= -9..9, transparency= .1),      
    plots:- pointplot3d(
        [[6, -8, f(6, -8)]], color= red, symbol= solidcircle,
        symbolsize = 18
    ), 
    plots:-spacecurve(
       [cos(t), sin(t), 1 - 12*cos(t) + 15*sin(t)], t = 0 .. 2*Pi,
       color= red, thickness= 3
    ),
    orientation= [-15, 68, 5],
    view= [-4.2..8.2, -8.2..4.2, -100 .. 100]
);

If you remove the local in the definition local ode_type:= module(), then it works, and you get the error message that you should get.