The expand command with split the integrals. You could do this:

Int(exp(-I*(A*x+B*y)), [x= -a...a,  y= -a..a]):
convert(value(expand(%)), trig);
 

I think that the HSV color scale provides a more-intuitive way than RGB of translating reals from bounded intervals to colors. The below code is a variation of Acer's, but I scale along H from 0 to .85 and keep S and V constant at 1. H values higher than about .85 look closer to red than violet to me. 

I have also provided what I think is a more intuitive way to use value to get the numeric values.
 

Download TypesetHSV.mw

It can be done like this:

rangetolist:= (x::range(realcons))-> [seq(x)]:

If you want to interpret the instructions literally, replace range(realcons) with (integer..integer).

It might be possible if the expressions are polynomials. See ?solve,parametric. The algorithm used by this command is extremely complicated, and results often take a long time to compute.

Using RealDomain seems to interfere with the ability of simplify to process assumptions. So, don't use it. I don't know exactly why it interferes, but I do know that RealDomain alters the definitions of both sqrt and simplify.

AFAIK, testeq is a probabalistic test that uses modular arithmetic. Because of the modularity, I don't think that it can ever process assumptions. You should replace testeq with is, which is a deterministic test that uses ordinary complex-number arithmetic and is the main command for deriving true/false conclusions from assumptions.

In the example below, S is the summand, which I've represented as an m-ary function:

n_0:= 4:  m:= 3:
add(S(seq(i +~ 1)), i= Iterator:-Combination(n_0, m));
      S(1, 2, 3) + S(1, 2, 4) + S(1, 3, 4) + S(2, 3, 4)

The inner seq will always sequence the combinations in order.

What you referred to as "all monotonic lidyd [lists] of size m of [1..n_0]" is usually called all m-combinations of {1..n_0}.

An older-style command with the same result is

add(S(i[]), i= combinat:-choose(n_0, m));

The function osum could be written like

osum:= (S, k::(name= [nonnegint, nonnegint]))->
    local j;
    add(eval(S, lhs(k)= j+~1), j= Iterator:-Combination(rhs(k)[]))
:

And called like

osum(i[1]^2+i[2]^2+i[3]^2, i= [4,3]);
 

You could do this:

j:= piecewise(x < 0, x^2 + 1, 0 <= x, x - 1);
plot(j, x= -3..3, discont);

If I do as above, I can right click on the output of the first line, then "Open Context Panel", then "Plot Builder".

The error message is fairly clear: In order to get a numeric solution for x, you need to supply numeric values for h, kpT[A], omega[A]. This is also mathematically clear. It seems unlikely that you're looking for a symbolic solution for x, because you've specified such a narrow numeric range for it.

I don't know what version of Maple you're using, but on mine it is already on the plot context menu. Right click on a plot for the context menu, then Symbol (7th item from the top in my version), then Symbolsize (bottom of the submenu), then enter the number.

The following procedure returns a root-bounding interval on which the original function is continuous, and it doesn't depend on the function or its numerator being polynomial.

p:= (x^3 - 13*x^2 + 55*x - 73)/(x - 1)
:
RootBoundwithDiscontinuity:= proc(f::algebraic, ab::name=range(realcons))
local A,k;
    for A in indets(plot(f, ab, discont), Matrix) do
        for k to upperbound(A,1)-1 do
            if A[k,2]*A[k+1,2] <= 0 then return [A[k,1],A[k+1,1]] fi
        od
    od;
    return
end proc
:
RootBoundwithDiscontinuity(p, x= -9..9); 
     [2.60815727428114, 2.64908352866142]

In addition to what Acer showed you, it's possible to combine any number of points, lines, and curves in a single plot command. Here's one with numerous bells & whistles.

L:= y = -3/2*x + 2; Pt:= [2,1]:
m:= -1/implicitdiff(L, y, x);
Lperp:= y = m*(x-Pt[1]) + Pt[2];
PtInt:= eval([x,y], solve({L, Lperp}));
plot(
    [rhs(L), rhs(Lperp), [PtInt], [Pt]], x= -5..5,
    style= [line$2, point$2], 
    symbolsize= 18, symbol= [solidcircle, soliddiamond],
    scaling= constrained, color= [red, "DarkGreen", black, red],
    view= [-5..5, -5..5], thickness= 0, size= [900$2],
    legend= [
        typeset(L, `\t`), 
        typeset(Lperp, ` `, `#mo(⊥)`,`\t`),
        typeset(PtInt, `\t`),
        Pt
    ],
    labels= [x,y], 
    caption= "\nThe given information is in red.", 
    captionfont= [Times, 18],
    title= "Perpendicular Lines\n", titlefont= [Helvetica, Bold]
);

Here's a different way that uses neither complex arithmetic (like VV's) nor an idiosyncratic fractal-specific syntax (like Tom's). I have nothing against those methods; I just think that this real-arithmetic-only approach provides another valuable way of understanding what's happening.

N:= 6: #number of iterations
#Rotate point P by angle theta around point Q:
Rot:= theta-> (P,Q)-> 
    local R:= P-Q, s:= sin(theta), c:= cos(theta): 
    <R[1]*c - R[2]*s | R[1]*s + R[2]*c> + Q
:
Rot120:= Rot(2*Pi/3): #the only rotation that we need
#Select a point on line segment PQ:
BetweenPt:= (P,Q,r)-> P*(1-r) + Q*r:
#Replace line segment __ with _/\, i.e., put a notch in the middle:
Notch:= (P,Q)-> 
    local P1:= BetweenPt(P,Q,1/3); 
    (P, P1, Rot120(P,P1), BetweenPt(P,Q,2/3))
: 
Pts[1]:= <seq(<cos(2*Pi/3*k) | sin(2*Pi/3*k)>, k= 0..3)>:
for k to N do
    Pts[k+1]:= 
        <seq(Notch(Pts[k][j], Pts[k][j+1]), j= 1..3*4^(k-1)), <1|0>>
od:
plots:-display(
   seq(plot(Pts[k])$4, k= 1..N+1), insequence, 
   thickness= 0, scaling= constrained, axes= none
);

I made 4 changes to your worksheet. The change to the plot command is syntactically required.  The other 3 changes are just to increase the speed of numeric integration:

1. Digits:= 15: #Increase from 10
2. E:= (x,t)-> Int(  #Use uppercase I.
       exp((-esp*w^4+Disp*w^2+k)*t)*cos(w*(x+v*t))/Pi,
       w= 0..infinity,
       epsilon= 1e-7  #Added for speed.
   );
3. #Same changes with larger epsilon:
   u:= (x,t)-> Int(E(x-xi, t)*f(xi), xi= 0..2e4, epsilon= 1e-5);
4. plot(u(1500, t), t= 0..6e4, numpoints= 100);

This is the plot I got for numpoints= 10:

You'll need to wait several minutes for numpoints= 100.

In the worksheet below, the procedure F represents your 2nd-order ODE as a 1st-order system in functional form. The procedure Heun specifies the iteration for what I believe your instructor means by "modified Euler". However, if this formula seems unfamiliar to you, then specify the formula that you're supposed to use, and this worksheet can be adapted by changing that 1 line.
 

Download PredictorCorrector.mw

If the integration is in rectangular coordinates, as it is in the example that you present, then a plot of the region of integration can be automatically generated from the inert form of the integral, like this:

J:= VectorCalculus:-int(x*y, [x,y]= Triangle(<0,0>, <1,0>, <0,1>), inert):
plot3d(0, op([1,2], J), op(2,J), orientation= [180,0,180]);

Then the integral can be computed by

value(J);