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);

There's no need to create a 60,000 point spline! Just plot the points:

X:= ExcelTools:-Import("E:/ct.xls", "Sheet1"):
Y:= ExcelTools:-Import("E:/ct.xls", "Sheet2"):
plot(
    <Y | X[..,2]>, 
    labels= ['t', 'concentration'], labeldirections=[horizontal, vertical]
);

Both parts a) and c) can be answered with the same vector cross product. This can be accessed by loading the LinearAlgebra package and using the &x operator. It's as simple as this:

restart:
with(LinearAlgebra):
V:= <x, y, f(x,y)>: P:= <5,10,0>: Q:= <-1,7,1>: R:= <2,1,4>:
#Order doesn't matter as long as all points are used:
solve((nv:= (P-Q) &x (R-Q)).(P-V), {f(x,y)})[]; #linear function
Norm(nv,2)/2; #area of triangle

The cross product of the direction vectors formed from any two distinct pairs of the triangle's vertices is a normal vector (nv) to the plane (i.e., graph of the linear function) containing the triangle. Thus, the plane consists of all direction vectors emanating from any vertex of the triangle (I chose P) whose dot (or inner) product with nv is 0. The 2-Norm (or magnitude or length or Euclidean norm) of nv is twice the area of the triangle.

I think that the problem occurs because Maple is confused by possible situations where if certain parameters were equal then the geometric multiplicity of certain eigenvalues would change. A workaround for this is to change the parameters to "safe", "generic", and "independent" transcendental constants, compute the eigenvectors, then change the constants back to the parameters. Like this:

restart:

phi:= Pi:
A := <
          <0, 0, exp(I*k1) + m1, exp(I*k2) + m2>|
          <0, 0, exp(I*phi)*(exp(-I*k2) + m2), exp(-I*k1) + m1>|
          <exp(-I*k1) + m1, exp(-I*phi)*(m2 + exp(I*k2)), 0, 0>|
          <exp(-I*k2) + m2, exp(I*k1) + m1, 0, 0>
     >:
Subs:= [k1= exp(3), k2= exp(sqrt(2)), m1= sin(1), m2= cos(sqrt(2))]:
EV:= [LinearAlgebra:-Eigenvectors(eval(A, Subs))]:
eval(EV, (rhs=lhs)~(Subs)); simplify(%);

Doing this gives you two eigenvalues, each associated with two linearly independent eigenvectors. A quick visual scan suggests that these are the same as reported by Mathematica, although it may be a challenge to prove that conclusively. A starting point may be to have Mathematica convert its results to sin/cos form.

I substituted for all four parameters, but I suspect that it would still work with substituting for fewer.

Something pretty close to that can be done like this:

plot3d(
    1e-2*abs(sin(6*t)*cos(6*x))*exp(-x-t), x= 0..1, t= 0..1,
    axis[1,2]= [
        tickmarks= [
            0= "0", 1="1", seq(.1*k= sprintf("%3.1f", .1*k), k= 1..9)
        ]
    ],
    axis[3]= [
        tickmarks= [
            seq(
                1e-3*k= typeset(sprintf("%3.2f x ", k), `10`^`-3`),
                k= 1..8
            ),
            seq(1e-3*(k+0.5) = "", k= 1..8) #subticks
        ]
    ],
    axesfont= [Helvetica, bold, 9],
    labels= [` x`, `t    `, `z   `], labelfont= [Helvetica, 32]
);

The "square" in your title suggests that you may only be interested in solutions that make use of repeated squaring. Acer's last paragraph is about Maple's built-in modular repeated squaring command: 

A &^ B mod 1000

If you're interested in detailed studying of repeated squaring from first principles, then here's some code to peruse:

restart
:
#Computes A*B mod b^d by convolution of lists of digits.
`&*`:= proc(A, B, b:= 10, d:= 3)
local c:= 0, i, k; #c is "carry" to the next digit
    [seq(irem(c + add(A[i]*B[k+1-i], i= 1..k), b, 'c'), k= 1..d)]  
end proc
:       
#Computes A^B mod b^d by repeated squaring.
`&**`:= proc(A::posint, B::posint, b:= 10, d:= 3)
local 
    sq:= convert(b^d+irem(A, b^d), ':-base', b)[1..d],
    r:= [1, 0$(d-1)], bits:= B, `&*`:= (A,B)-> :-`&*`(A, B, b, d)
;
    while bits <> 1 do 
        if irem(bits, 2, 'bits')=1 then r:= r &* sq fi;
        sq:= sq &* sq
    od;
    add(r&*sq*~b^~($0..d-1)) #Recompose digit list to integer output.
end proc
:

To use it for the problem at hand:

245 &** 272

As far as I know, that can't be changed. Even overloading > automatically overloads <, and likewise for >= and <=. Even unevaluation quotes won't stop the switch. For example, try ' ' a > b ' ' (that's two pairs of single quotes, not one pair of double quotes).

I'll treat here only the k=1 case. Other specific positive integer values of k can be handled similarly.

Download sincossum.mw

A good command for plotting the curve of intersection between two surfaces is intersectplot. Below, I apply it pairwise to the set of three planes.

plots:-display(
    plots:-implicitplot3d(
        [sys[]], (x,y,p)=~ 0..10, 
        color= COLOR~(RGB, [.5,0,0], [0,.5,0], [0,0,0.5]),
        grid= [3$3], transparency= 0.2, thickness= 0),
    seq(
        plots:-intersectplot(
            P[], (x,y,p)=~ 0..10, thickness= 9, color= magenta
        ),
        P= combinat:-choose(sys, 2)  #pairwise
    ),
    orientation= [160, 70], projection= .8, 
    lightmodel= light4, axes= frame, tickmarks= [3$3]
);