If you literally want to plot just the points, it can be done by

plots:-pointplot3d(
   eval([seq(seq([a+i*h, c+j*k, U[i,j]], i= 1..nX-1), j= 1..nY-1)], soln),
   symbolsize= 16, labels= [x,y,U]
);

But I think that you may be more interested in plotting the surface that can be interpolated from the points. That can be done by

plots:-surfdata(
   eval(Matrix((nX-1, nY-1), symbol= U), soln),
   a+h..a+(nX-1)*h, c+k..c+(nY-1)*k, labels= [x,y,U]
);

Neither is done literally with the plot3d command. Rather, closely related commands are used.

Simply define the following procedure:

`print/Unit`:= u-> [[u]]:

You do not need to ever explicitly use this procedure; merely having defined it is enough. This only affects the prettyprinted display of the units. It has no effect on computation.

p_yx:= diff(p(x,y), y, x):
thaw(value(subs(p_yx= freeze(p_yx), %eliminate({Eq1, Eq2}, p_yx))));

Maple's mod package (see ?mod) has everything that you need for working with polynomials over finite fields. Vv's Answer seems to imply that mod only works with finite fields of prime order. That's not true. Here's an example of computing the gcd of two polynomials over GF(3,2):
 

Download GF.mw

This quick Answer is all that I can manage right now: A parametrized RootOf is intended for expressing inverse functions symbolically. If you're going to evalf it for numerics, then just use fsolve instead. A procedure that uses fsolve can be numerically integrated; indeed, it is common to do so.

There are many ways to do this. The way that I prefer is to use a single plot command with a seq-loop inside instead of a plot command inside a for-loop. Then no display is needed.

HA:= [1,2,3]:
n:= nops(HA):
plot(
   [seq(eval([f[1], f[2], f[3]], R= HA[i])[], i= 1..n)], 
   x = 0 .. 1, size = [300, 250], labels = [y, T(y, t)], 
   titlefont = ["Arial", 10, Bold], 
   labelfont = ["Arial", 10, Bold], 
   labeldirections = [horizontal, vertical], 
   axes = boxed, color= [blue $ n, green $ n, red $ n], 
   linestyle = op~([[solid, dash, dot] $ n]), thickness = 3
); 

You may want to include a legend, if you can come up with some succinct names for each of the nine curves.

The easiest way is to do this in a single execution group:

st:= time():
# Code that you want to test...
# ....
time() - st;

In other words, you mark the starting CPU time and subtract it from the finishing CPU time.

If your code can be expressed as the evaluation of a single expression---which can always be achieved by putting the code into a procedure---then you can use the simpler

CodeTools:-Usage(expression);

This command returns the value of the expression and prints timing and memory usage data as a side effect.

From some preliminary plot experimentation, I concluded that for real x and y, the eigenvalues are all real in the interval [0,7]. Then an Explore command lets you plot the eigenvalues while using sliders to control the values of x and y. I included code so that imaginary eigenvalues or those outside the given real range would, if any are found, would be reported with print statements.

CP:= (x,y)-> 
   -(lambda-2)^2*(544-lambda^9-328*lambda^7-16*cos(x-y)+6*cos(x-y)*lambda+
   6*cos(x+y)*lambda-408*cos(y)*lambda^2-408*cos(x)*lambda^2+
   540*cos(y)*lambda+540*cos(x)*lambda-28*cos(y)*lambda^4-28*cos(x)*lambda^4+
   152*cos(y)*lambda^3+152*cos(x)*lambda^3+2*cos(y)*lambda^5+2*cos(x)*lambda^5-
   256*cos(y)+17232*lambda^2-256*cos(x)-7852*lambda^5+17768*lambda^4+
   2088*lambda^6+28*lambda^8-16*cos(x+y)-23632*lambda^3-5844*lambda)
:
EV:= proc(x,y)
local r;
   if evalf([x,y])::'list(numeric)' then
      r:= [fsolve(CP(x,y), 'complex')];
      if not r::['realcons' $ 11] or min(r) < -.1 or max(r) > 7 then
         print("Deviant eigenvalue detected:", r, "at", [x,y]) 
      end if;
      r
   else
      [%fsolve(CP(x,y), 'complex')]
   end if
end proc:

Explore(
   plot(
      `[]`~(EV(x,y), 0), 
      style= point, symbol= diamond, symbolsize= 16, 
      axes= box, view= [-1..7, -.1..0.1], size= [1500,200]
   ), parameters= [x= -Pi..Pi, y= -Pi..Pi]
);

The root of your problem is the failure to use unapply and really has nothing to do with piecewise, as detailed in Kitonum's Answer. I just want to point out a simpler way to write piecewise expressions: Any piecewise of the form

piecewise(c, e1, not c, e2)

can be expressed as

piecewise(c, e1, e2).

And more generally, any piecewise of the form

piecewise(c[1], e1, not c[1] and c[2], e2, ....,  not `or`(c[1], ..., c[k-1]) and c[k], ek, ...)

can be expressed as

piecewise(c[1], e1, c[2], e2, ..., c[k], ek, ...).

In other words, they are processed left to right such that every condition already includes the negation of all conditions to the left of it.

You have

Sol[1]:= [fsolve([...] union {...})];

That needs to be

Sol[1]:= [fsolve({...} union {...})];

You need to pay attention to whether the interval of integration crosses vertical asymptotes and restrict x accordingly:

int((s^3-5*s^2+1)/(s+2)^2, s = 0 .. x) assuming x > -2;

Have a look at ?Statistics,KernelDensity.

The instructions don't say to replace solve with evalf , they just say to use evalf. The intention was that you first get exact expressions for the critical points with solve. This involves solving some quadratic equations, so there are some square roots in these answers if you input the equations with exact fraction coefficients rather than decimals. Then you can get some decimals for these with evalf. Try this

EquationSet:= {
   diff(x(t),t) = -y(t) - z(t), 
   diff(y(t),t) = x(t) + y(t)/5, 
   diff(z(t),t) = 1/5 + (x(t) - 5/2)*z(t)
};
Solut:= [solve(rhs~(EquationSet), {x,y,z}(t), explicit)];
evalf(%);

Your matrix essentially specifies a function from R^2 to R^4. So there are 6 dimensions. One possible way to represent this is as four surfaces---one for each coordinate function. These can be shown as a single 3D plot of four surfaces or as four separate frames arranged as a 2 x 2 array. No matter how you represent it as a plot, you'll need ranges for x and y.

Somewhere in your input expression---likely in several places---you have used square brackets [.] were you should've used round parentheses (.). All algebraic grouping in Maple, no matter how deeply nested, is done with parentheses.