You can't use square brackets as parentheses---as an outer level of grouping, for example. You can only use round parentheses for algebraic grouping. So, your expression should be like this:

solve({((alpha[1]-alpha[2]*lambda)*sqrt(x)+p[2]*lambda)*(k[1]*(1-lambda^2)+2) = p[1]*(2*k[1]*(1-lambda^2)+2), ((alpha[2]-alpha[1]*lambda)*sqrt(x)+p[1]*lambda)*(k[2]*(1-lambda^2)+2) = p[2]*(2*k[2]*(1- lambda^2)+2)}, [p[1], p[2]]);

P:= Matrix([''[1,4],[1,5],[2,6],[3,7]'' $ 3]):
x:= 3*Statistics:-MovingAverage(P[..,1], 3):
X:= 3*Statistics:-MovingAverage(P[..,1]^~2, 3):
y:= 3*Statistics:-MovingAverage(P[..,2], 3):
Y:= 3*Statistics:-MovingAverage(P[..,2]^~2, 3):

select(has, f, u);

select(has, f, x);

Use this procedure:

MatrixToStrings:= proc(M::Matrix(integer))
uses LA= LinearAlgebra;
local k,x;
     [seq(cat(seq(sprintf("%d",x), x= M[k,..])), k= 1..LA:-RowDimension(M))]
end proc:

Use add, not sum.

add(add(add(a_(i,j,k)*d(i, j, k), k = 0 .. 1), j = 0 .. 1), i = 0 .. 1);

The reason that it doesn't work with sum is that i, j, and k have the value 1 from some previous computation. Using add, it will ignore previous values of the index variables. But sum---like the vast majority of Maple procedures---evaluates its arguments before they are passed.

f:= n-> a^n + b^n + c^n:
Sol:= {solve({f(1)=1, f(2)=2, f(3)=3}, {a,b,c})}:
nops(Sol);
                               1
evala(eval(f(5), Sol[1]));
                               6

Here's a solution technique that tends to give a large number of zeros in the solutions vector, and a large number of integer entries.

I will assume that your 500x1300 matrix contains floating-point data. If it contains integer data, then you can skip the Integerize step below.

The first is convert the problem to an integer problem by multiplying each row by its least common denominator. The procedure Integerize does that.

Download Integerize.mw

This is probably not more efficient or direct, but it does look neater to me and easier to understand:

tmp:= seq(phi[j-1]*~[1-p[j], p[j]], j= 5..2, -1);
Vector(op~([tmp]));

A *~ B

See ?elementwise

It takes several short steps to remove all of Maple's knowledge about gamma. In Maple 16, do these:

restart:
constants:= ({constants} minus {gamma})[]:
`evalf/gamma`:= proc() end proc:
`evalf/constant/gamma`:= proc() end proc:
unprotect(gamma);

Note that this will affect certain special functions, such as Psi, that depend on gamma. The effect will ensue even when those functions are being used (unbeknowst to the user) in the background, rather than being explicitly invoked by the user.

Consider the following plot:

plot([x, 1/3 - (1/3)*x], x= 0..1);

Since the two lines intersect between x=0 and x=1, there is no meaningful way to order the expressions; the is will return false no matter which direction you put the inequality.

sys:= {diff(y(t),t) = 2 - t^2/2 - 1/4 * Y(t), diff(Y(t),t)= y(t), Y(0)=0}:
eval(y(t), dsolve(sys, {Y(t), y(t)}));

Here's a crude procedure for it:

Download PowerMethod.mw

You can do it with LinearAlgebra:-QRDecomposition or LinearAlgebra:-LinearSolve(..., method= QR).

There is some bug here; I don't know specifically where. The bug is manifested regardless of whether there are assumptions.

Here's a workaround: convert to a list and then back to a Vector.

k:= Vector[column]([W, (1/3)*W, (1/3)*W, (1/3)*W]);
Vector(sort(convert(k, list), (a,b)->is(a>b)));