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

Yes, you can substitute the identity function or any other function:

eval(Sol, _F1= (x-> x));

where Sol is whatever you named the solution.

Your curve can be parameterized to a space curve, which can be plotted with tubeplot, which unlike spacecurve does accept the colorscheme option with suboption xyzcoloring.

restart:
S1:= x*y+y*z+x*z = 0:
S2:= x+y+z = 5:
E:= eliminate({S1,S2}, z);
Polar:= {x= r*cos(t), y= r*sin(t)}:
R:= solve(eval(E[2], Polar), r);
XYZ:= combine(eval[recurse]([x,y,z], Polar union E[1] union R[2])); 
plots:-tubeplot(
    XYZ, t= -Pi..Pi, radius= .2, numpoints= 300, 
    colorscheme= ["xyzcoloring", `*`], style= surface, scaling= constrained
);