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

There are two ways: You can make a new vector of the results, or you can store the modified elements back in the original vector. This latter form is called inplace operation. Understanding the difference between inplace and non-inplace operation is often important when working with rtables (Vectors, Matrices, and Arrays).

The selection operation is

Lambda[2..4] =~ 0

This makes a new vector, which I've not named. You can assign it to a name, convert it to a set or list, etc. 

MyEqnSet:= {seq(Lambda[2..4] =~ 0)}  #a set of equations

The inplace operation is

Lambda[2..4]:= Lambda[2..4] =~ 0;

To specify more than one PDE, put them into a set or a list. A set is enclosed with curly braces:

MySys:= { PDE1, PDE2 };

A list is enclosed with square brackets:

MySys:= [ PDE1, PDE2 ];

The elements of a list or set must be separated by commas.

Most high-level commands, including those shown in your worksheet, will accept either form. In the case of a set, Maple will rearrange the elements into Maple's own fixed order.

plots:-fieldplot3d([-x, y, cos(z)], x= -Pi..Pi, y= -Pi..Pi, z= -Pi..Pi);

In addition to all the usual options for 3d plots, there are a vast number of options to control the length, thickness, color, style, placement, and number of arrows (see help page ?fieldplot3d).

Add the clause assuming y > 0 to your simplify command.

If you read the algorithm-details page that your linked page links to, you'll see that they start from a massive database of precomputed values. That database likely takes up most of a large disk drive. If you want a good estimate that'll be returned instantly:

Digits:= 15:
trunc(evalf(Li(6469693230)));

I know virtually nothing about integral equations, but the following totally ad hoc seat-of-the-pants method converges quickly to a solution:
 

Download IntEq.mw

Steps:

1. Make g simply the piecewise expression, without the x-> part.
2. Change the range to x= -5..5.
3. Change discont= true to discont= [showremovable].
4. Add undefined as a last argument to piecewise. (This is so that it doesn't default to 0 for the unspecified gaps.)

Either replace the keyword trig with arctrig, or remove it entirely.

I can find no help page for simplify(..., arctrig), but showstat(`simplify/arctrig`) shows that it exists.

There's a simple closed-form solution for arbitrary a, b[0], b[1], b[2], and u[1]. Take the expanded u[10], for example, and collect wrt the bs. The leading term is a^9*u[1]. The coefficient of each of the bs is a geometric sum of powers of a^3.

I'm recovering from the flu at the moment, and I don't have the energy to work out all the details at the moment.

The whole group can be created with a single cat command; no loop is needed:

ourfile:= cat("C:/Users/Tamour Zubair/Desktop/AA/T[", 1..5, "].dat");