You're simply not correctly quoting the file name or including the file extension. Try this:

#Create list L and Matrix A before doing the save.

#Verify L is a list and A is a Matrix:
type([L,A], [list,Matrix]); #Should return true
save L, A, "mysavefile.mpl";
restart;
read "mysavefile.mpl": #End with a colon to avoid lengthy output.
type([L,A], [list,Matrix]);

If the last command returns true then it worked; if it returns false then report back here and we'll fix it.

I think that the above addresses all your questions, but to be explicit:

• Does the file need to be created first, before being saved to?

No, and if it has been created, the save command will simply overwrite it.

• Where should that file be placed? Desktop? Folder?

It makes no difference, as long as you save to and read from the same place. In the above example, I simply used whatever the default folder was.

• How do I read from the file after I have saved?

Use the read command.

• Can I save both L and A to the same file?

Yes.

• Do I need 2 commands?

No. If you want them in the same file, then you must save them with a single command.

I wrote a module to do it:

EnvelopePlot:= module()
export
   Extract:= (P::specfunc(PLOT))->
      [plottools:-getdata(P, "curve")][1][-1],
   
   Upper:= proc(C::Matrix)
   local X:= C[..,1], Y:= C[..,2];
      if not LinearAlgebra:-Equal(X, sort(X)) then
         error "CURVE isn't sorted wrt its 1st coordinate."
      end if;
      <X |
         CurveFitting:-ArrayInterpolation(
            C[select(
                 k-> Y[k-1] <= Y[k] and Y[k+1] <= Y[k],
                 [$2..op(2, [rtable_dims(Y)][1])-1]
              ), ..],
            X,
            method= spline
         )
      >
   end proc,

   Lower:= (C::Matrix)-> 
      < C[..,1] | -Upper(< C[..,1] | -C[..,2] >)[.., 2] >
;
end module:   

And here's an example of its use on an odeplot:

S:= dsolve({diff(f(x),x) = diff(cos(16*x)/(1+x^2), x), f(0)=1}, numeric);
P:= plots:-odeplot(
   S, [x,f(x)], x= -2*Pi..2*Pi, 
   numpoints= 2^10, color= black, legend= ["Signal"]
):
C:= EnvelopePlot:-Extract(P):
plots:-display(
   P, 
   plot(EnvelopePlot:-Upper(C), color= red, thickness= 3, legend= ["Upper envelope"]),
   plot(EnvelopePlot:-Lower(C), color= blue, thickness= 3, legend= ["Lower envelope"],
   axes= boxed)
);

The maximum number of objective-function evaluations is set by default to 50. You can change that by passing option evaluationlimit= 99 (or whatever number) to Maximize. I did this for your PD(10, 10, 10, 200), and there was no problem; the max number of evaluations that was actually used was 53.

Regarding the gap problem: I suspect that there's a bug in your code, but I'm not sure.

_U1 is a "bound", yet global, variable. It is the variable of integration. When Maple needs to make up a variable name to place in output, the name is of the form underscore-letter-integer. You can use subs to replace the name with another if you want. Note that t is used in the integral, so replacing _U1 with t would lead to an incorrect result. 

Z:= sqrt(1 - r^2*sin(theta)^2):
plots:-display(
   plot3d([r, theta, Z], r= 0..1, theta= -Pi..Pi, coords= cylindrical),
   plot3d([r, theta, -Z], r= 0..1, theta= -Pi..Pi, coords= cylindrical),
   plot3d(
      [1, theta, z], theta= -Pi..Pi, z= eval(-Z..Z, r= 1), coords= cylindrical, 
      style= wireframe, color= purple, thickness= 4
   ), 
   scaling= constrained, labels= [x,y,z]
);

#Volume integral in cylindrical coordinates:
Int(r, [z= -Z..Z, theta= -Pi..Pi, r= 0..1]):
% = value(%);

#Volume integral in cartesian coordinates:
Int(1, [z= -sqrt(1-y^2)..sqrt(1-y^2), y= -sqrt(1-x^2)..sqrt(1-x^2), x= -1..1]):
% = value(%);

It's clear by inspection that (x,y) = (0,0) is a solution for any t. It can be difficult for numerical methods to avoid converging on a trivial solution. Do you expect other solutions, and if so, can you give ranges in which you expect them to fall? They'll need to be ranges that don't include 0.

If we use a definite integral with a variable upper limit and some reasonable assumptions, then we'll get the obvious result based on the binomial expansion of (x-a)^m.

restart:
value(int((x-a)^m/x, x= a..x)) assuming x > a, a > 0, m::posint;
#Terms without x can be discarded as constants of integration:
select(has, %, x);
#Convert GAMMA to factorial:
convert(%, factorial);
subs(_k1= k, %);
simplify(%);

-x^m*(sum((-1)^k*a^k*x^(-k)/((-m+k)*k!*(m-k)!), k = 0 .. m-1))*m!+a^m*(-1)^m*ln(x)

What you are calling the "gradient" of a vector field is called the Jacobian in English. You simply need to change Gradient to Jacobian.

In English, as far as I know, one only computes the gradient of functions from R^n to R^1. The Jacobian is a generalization of the gradient to functions from R^n to R^m. The distinction is only linguistic and a matter of convention; there's no practical difference between first derivatives, gradients, and Jacobians.

When I used what seemed to me to be more-realistic values for the constants, Maple wasn't as compliant as it was for Kitonum's values of the constants.

Numeric integration over a hyper-rectangle is more robust than integration over regions with variable inner limits of integration. I can see immediately that your region of integration is a cylinder of height t over the unit disk. So a conversion to cylindrical coordinates will make the region of integration a hyper-rectangle.

Download NumIntCylCoord.mw

I can confirm that value to 9 digits by using the original coordinates, performing the innermost integral symbollically, and then the two outer integrals numerically. I don't think that this method is as robust numerically as the conversion to cylindrical coordinates.

A conditional solution will either be a piecewise branch with an inequality condition (if you designated which names were to be considered parameters), or it'll contain inequalities (if you didn't). The presented case is the latter. To separate the solutions which contain inequalities from those that don't:

S:= solve({2*x+y+z = 3, x^2+y^2+z^2 = 3}, [x, y, z], real, parametric);
selectremove(x-> membertype({`<`, `<=`}, x), S);

To just get the unconditional solutions, change selectremove to remove.

For what it's worth:

I changed your initial condition so that it didn't reuse the independent variable x. So, I made it u(a, 1/a) = 0. This eliminated the error message, so I guess that that is the preferred form for initial conditions; however, I let the pdsolve command run 11 hours without getting any result, then I killed it. Like a zombie, it started running again. I eventually needed to kill my entire Maple session. This is in Maple 2016.

If you include a trivial inequality in the system, solve seems spurred to return only real solutions without the need to specify any other options, load any packages, or use any assumptions. This is even true if the inequality has no relation to the other equations in the system, as in

solve({x^2+y^2+z^2=3, x+y+z=3, w > 0}, {x,y,z});

     {x = 1, y = 1, z = 1}

I obtained the above result both with and without RealDomain.

Here's the procedure npar:

npar:= proc(p::procedure, $)
local L:= [op(1, eval(p))];
   nops(L) - `if`(` $` in L, 1, 0)
end proc;

This counts the parameters, not the arguments. What's the distinction? Parameters are an inherent property of a procedure; they're part of how it's written. Arguments are what's used when the procedure is called. The number of arguments to a given procedure may vary. If I declare

f:= x-> x^2;

then f has one parameter, x. If I call it using f(1,2), then it has two arguments for this call. This will not generate an error; by default it's allowed to pass extra unused arguments to a procedure. If you want to insist that f be called with exactly one argument, then you need to write it like this:

f:= proc(x,$) x^2 end proc:

Now that $ will appear as ` $` in op(1, eval(f)), which is why my procedure npar checks for it and subtracts 1 if it's found.

Inside procedure lich, you use f as if it were itself a procedure. That's fine, but then you need to pass a procedure to f, not an expression. So, you just need to change the call to lich to

lich(x-> sin(x^2), 0, Pi, 10);

We usually don't take kindly to screenshots. In the future, please upload a worksheet using the green uparrow on the toolbar of the MaplePrimes editor.

To get longitude and latitude lines on a sphere, use plot3d with coords= spherical:

plot3d(
   [1, phi, theta], theta= 0..2*Pi, phi= -Pi/2..Pi/2,
   coords= spherical, scaling= constrained
);

You can control the number of lines (to some extent) with the grid= [m, n] option.