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

The Power method won't converge for this matrix because there are two eigenvalues with the maximal absolute value.

You should stop using the deprecated linalg package. Everything you need is in LinearAlgebra, for example Eigenvalues and Eigenvectors.

You just have some minor syntax problems:

1. You need an extra pair of square brackets to enclose the data.
2. You used two independent variables, x and v, in the curve.

CurveFitting:-LeastSquares(
     [
          [1.0, 2.33], [2.0, 0.626e-1], [3.0, -2.16],
          [4.0, -2.45], [5.0, -.357], [6.0, 2.21],
          [7.0, 2.75], [8.0, .636], [9.0, -2.45]
     ],
     x,
     curve = a+b*cos(x)+c*sin(x)+d*cos(2*x)+e*sin(2*x)
);

Note that this solution does not explicitly use the LinearAlgebra package. It can be modified to use LinearAlgebra:-LeastSquares if that's what you want.

You should stop using the deprecated linalg package in all your code.

Download Interpolate.mw

I used exact arithmetic in your worksheet, and I got all zeros at the end. To use exact arithmetic, I did three things:

1. I changed all numbers of the form 210e9 to 210*10^9
2. I changed all numbers of the form 0.5, 0.25, etc., to 1/2, 1/4, etc.
3. I removed all evalf.

Luckily, there was no significant increase in the computation time. With these changes, the setting of Digits is insignificant.

Here is the modified worksheet: Exact_arithmetic.mws