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

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