Use one of the following:

plot3d(p1, -1 .. 1, -1 .. 1);

expr:=piecewise(x^2+y^2 <= 1, x*y-y^2, 0);
plot3d(expr, x = -1 .. 1, y = -1 .. 1);

plot3d('p1(x,y)', x = -1 .. 1, y = -1 .. 1);

The command to plug n=1 into the expression S is:
eval(S, n=1);

The command solve is to find the root(s) of the equation S=0, and the syntax solve(S, n=1) is wrong.
You should read (or at least browse) the Maple User Manual.

No need for assumptions, but instead:

- Don't use square brackets [ ... ],  they are list delimiters.
- Don't use inert Int , or use value(%) after that.
- (optional) use Maple (1D) input.

int( (x/L)^3 - 1, x=0..L );

           
                              

Yes, arrays of plots refuse to be exported (by context menu or by plottools:-exportplot).
Workaround:

restart;
p:=seq(plot(x^k,x=-1..1),k=1..4):
plotsetup(png,plotoutput="d:/temp/ex-array.png");
plots:-display(<p[1],p[2];p[3],p[4]>);
plotsetup(default);

subs(_Z9 = n, rhs(symbolic));
(of course the numeric suffix may differ).

 

When coords=cylindrical,  r  is ploted as a function of theta and z.
Even in your worksheet the notation is misleading; you use plot3d( f(r), theta=..., r=...)
but actually r should be z; it does not matter, of course, but your intentions are then not clear.

evalb(evalf[3](Pi=3.14));
                              true
evalb(evalf[4](Pi=3.14));
                             false

x0,y0,z0:=1,2,3:
eq:= u*~[cos(t),sin(t),0] + (1-u)*~[x0,y0,z0]:
plot3d( eq, u=0..1,t=0..2*Pi, scaling=constrained );

For the 2nd order derivative you must use

fdiff(dydata, [1,1], 0);

[Note that to compute this it would be much better to use the RHS of the ODE].

With your new constraints Maple cannot find a feasible integer point because such point does not exist. i.e. your constraints are incompatible.
You can see this using Maple (but also by hand):

The constraints are:

m := 9: n := 9: p := 3:
constr:= seq(add(Y[i,j],i=1..m)=1,j=1..n),
         add(X[i],i=1..m)=p,
         seq(seq( Y[i,j]<=X[j],i=1..m),j=1..n);

By symmetry we may take
X[1]:=1;X[2]:=1;X[3]:=1; for i from p+1 to m do X[i]:=0 od;

But Y[i,j]<=X[j] implies:

for i to m do for j from p+1 to n do Y[i,j]:=0  od: od;

and evaluating
constr;  # ==>

Y[1, 1]+Y[2, 1]+Y[3, 1]+Y[4, 1]+Y[5, 1]+Y[6, 1]+Y[7, 1]+Y[8, 1]+Y[9, 1] = 1, Y[1, 2]+Y[2, 2]+Y[3, 2]+Y[4, 2]+Y[5, 2]+Y[6, 2]+Y[7, 2]+Y[8, 2]+Y[9, 2] = 1, Y[1, 3]+Y[2, 3]+Y[3, 3]+Y[4, 3]+Y[5, 3]+Y[6, 3]+Y[7, 3]+Y[8, 3]+Y[9, 3] = 1, 0 = 1, 0 = 1, 0 = 1, 0 = 1, 0 = 1, 0 = 1, 3 = 3, Y[1, 1] <= 1, Y[2, 1] <= 1, Y[3, 1] <= 1, Y[4, 1] <= 1, Y[5, 1] <= 1, Y[6, 1] <= 1, Y[7, 1] <= 1, Y[8, 1] <= 1, Y[9, 1] <= 1, Y[1, 2] <= 1, Y[2, 2] <= 1, Y[3, 2] <= 1, Y[4, 2] <= 1, Y[5, 2] <= 1, Y[6, 2] <= 1, Y[7, 2] <= 1, Y[8, 2] <= 1, Y[9, 2] <= 1, Y[1, 3] <= 1, Y[2, 3] <= 1, Y[3, 3] <= 1, Y[4, 3] <= 1, Y[5, 3] <= 1, Y[6, 3] <= 1, Y[7, 3] <= 1, Y[8, 3] <= 1, Y[9, 3] <= 1, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0, 0 <= 0

As you can see, there are some "0=1" here !

You will have to reconsider the constraints.
 

Only a few corrections were needed.

p-Median-ok.mw

This should work.

plots[fieldplot3d]([2*x, 2*y, 2*z], x = -1 .. 1, y = -1 .. 1, z = -1 .. 1, grid = [5, 5, 5], arrows = SLIM,
colorscheme=["xyzcoloring", (x, y, z)->sin(x)*sin(y)*sin(z)]);

It is accepted syntactically but it seems to be ignored.
In the help:
"Most but not all plotting commands support this option".

They are not quite similar. The multi-ellipses are always convex while the lemniscates could be even nonconnected.
On the other side, when the radius tends to infinity, both become almost circles.

lemni:=proc(p,R:=1,M:=10)
local x,y, z:=indets(p,symbol)[1], P:=eval(simplify(p),z=x+I*y):
P:=evalc(P*conjugate(P)-R):
plots[implicitplot](P, x=-M..M,y=-M..M, numpoints=10000)
end:

lemni((z^7-1)/(z-1),1);
lemni((z^7-1)/(z-1),1000);

M:=LinearAlgebra:-JordanForm(A, output='Q');