A routine which would allow efficient iteration over an rtable (Array, Matrix, Vector) as well as allow simultaneous access to the indices as each element is accessed would be very nice. It's unfortunate that it does not exist; using side-effects with rtable_scanblock is inefficient and does not scale well.
You can map or Map over an rtable, but not get access to the indices as each element is dealt with. The same goes for syntax like,
for x in A do .... od: # for rtable A
As for the NULLs in your example, select() and its friends produce the same sorts of objects that they are passed. And it seems that the size of an Array is a quality preserved by select. So, some value is needed, and NULL is better than a default of 0, since 0 might have some special meaning for a given example of selection.
acer

You can write a routine to do this, using loops. Or you could have it use rtable_scanblock, which can be made to walk over Arrays, Matrices, and Vectors in all sorts of useful ways. One neat aspect is that, while walking such a structure, rtable_scanblock has access to the indices and so can save them elsewhere as a side-effect.
For example,
FindIndex := proc(comparison,A)
local checkindex,G;
checkindex := proc(val,ind)
if comparison(val) then
G[val] := op(ind);
end if;
NULL;
end proc:
rtable_scanblock( A, [rtable_dims(A)], checkindex, [[1],A[1]]);
convert(G,list);
end proc:
foo := Array( 1..10 , i -> i ) * 2;
FindIndex( x -> evalb( x mod 3 = 0 ), foo );
You could make it more general, depending on whether you want it to handle rtables of several dimensions or to return something other than a list.
acer

Consider this sort of looping,
for i from 1 to 100 do
cat( "p", convert(5000+i,string) );
od;
So, you can produce strings, for names of files, programmatically.
acer

It's not clear whether you expect to create elements which have 20 decimal places of information, or whether you just want the Matrix to be able to hold such numbers.
Matrix(2,2,(i,j)->evalf[20](1/(i+j+1)),datatype=sfloat);
If the datatype of the Matrix is 'sfloat', which stands for software float, then the entries can have any number of digits. The precision of what can be assigned into the entries, after creation of such a Matrix, would depend on the environment variable Digits.
See the help-page, ?UseHardwareFloats .
Other ways to get a datatype=sfloat Matrix are,
restart:
UseHardwareFloats:=false:
Matrix(2,2,(i,j)->evalf[20](1/(i+j+1)),datatype=float);
MatrixOptions(%,'datatype');
restart:
Digits := 20:
Matrix(2,2,(i,j)->evalf(1/(i+j+1)),datatype=float);
MatrixOptions(%,'datatype');
acer

Sorry, I don't know why, but I missed that L is L(x) a function of vector x.
Let's have another shot at it, then.
There seems to be a missing multiplication sign, in your definition of L, in the [1,1] entry. Should there be a * between exp((x[1]+x[2])/K) and (1-1/a) ?
L:= Matrix([[s[1]*exp((x[1]+x[2])/K)*(1-1/a),F],[s[1]*exp((x[1]+x[2])/K)*1/a,s[2]]]);
VV := Vector([x[1],x[2]]);
Z := L.VV - VV;
e1 := isolate(Z[2],x[1]);
solve(eval(Z[1],e1),{x[2]});
sol2 := op([%][2]); # pick the non-trivial solution
sol1 := op(solve(eval(Z[2],sol2),{x[1]}));
op(solve(Z[2],{x[1]}));
newZ := eval(Z,{sol1,sol2});
# Now, how to simplify that to zero?
# Are there some assumptions on the other unknowns, to give
# simplify() a helping hand?
# This is crude.
evalf[100]( eval( newZ, {K=1/2,s[1]=1,s[2]=2,a=3,F=1/2} ) );
acer

If L-IdentityMatrix(n) has full rank and if its determinant is not zero, then yes, only the trivial solution of the zero-Vector should get returned. I suspect that, in such a case, there is no other solution to be found, by any method.
But if L-IdentityMatrix(n) is not of full rank and has a determinant of zero, then NullSpace should give a basis for the solution space, or LinearSolve (or `solve` of your fully formed explicit system) should give parametrized solutions.
In other words, I don't think that the method is what's important here. What should be key is whether there exist non-trivial solutions (or not).
acer

Alternatively, without having to muck about with the mouse, or assistants,
MM := ImportMatrix("4399_p0000000.xls",datatype=anything);
Of course, you may need to make the string point directly to the .xls file with its full path, or to change currentdir() before making the call.
acer

Does this follow?
x = L . x
0 = L . x - Identity . x
(L-Identity) . x = 0
So, in Maple commands,
with(LinearAlgebra):
X := LinearSolve( L-IdentityMatrix(n), ZeroMatrix(n) );
where n is ColumnDimension(L). Or you could get a basis for the nullspace in which x lives, with,
NullSpace( L-IdentityMatrix(n) );
acer

This should work in Maple 10, as well as 11.
Open the ImportData assistant, either through the Standard GUI's menu bar (Tools -> Assistants -> Import Data), or using the explicit command ImportData() .
After doing that, I selected your .xls file. I changed the choice for the datatype from float[8] to 'anything', so that it could hold the first row of names.
Upon return from this assistant, I had a 2001x8 Matrix, which I assigned by issuing,
M := %;
acer

If you use t as the dummy variable of integration alongside t[0] or t[o], then t isn't really free. This is because names like t can refer to tables in Maple, and t[0] references the 0 index entry of t. Hence mixing t and t[...] is a muddle to be avoided, with semantic baggage that confuses Maple.
acer

On reflection, it may not be at all obvious to the new user how to store the data and the final iterates, so that they can be easily plotted.
restart:
N := proc(f)
local Xnew, X0, incX, A, B, k, df, count, Lreal, Lcplx, R;
R := Array(1..101*101,1..3);
df := diff(f, x);
# These two outermost loops are used to create the X0 initial points.
count := 0;
for A from 0 by 0.02 to 2 do
for B from -1 by 0.02 to 1 do
X0 := evalf(A + B*I);
# Now do ten Newton iterations, using X0
Xnew := X0;
for k to 10 do
incX := evalf(eval(f/df, x = Xnew));
Xnew := Xnew - incX;
end do;
count := count + 1;
R[count,1],R[count,2],R[count,3]:=evalf(A),evalf(B),Xnew;
end do;
end do;
# Return the results in lists suitable for pointplot3d.
Lreal := [seq([R[k,1],R[k,2],Re(R[k,3])],k=1..101*101)];
Lcplx := [seq([R[k,1],R[k,2],Im(R[k,3])],k=1..101*101)];
return Lreal,Lcplx;
end proc:
Lreal,Lcplx := N(sin(x)):
plots[pointplot3d](Lreal,axes=boxed);
plots[pointplot3d](Lcplx,axes=boxed);
acer

So, you had to use Newton's method to get from X0 to X1, X2, etc? If so, then what was the function? You must know, since you already solved question a). Someone else here also asked for this. No doubt it was laid out earlier, perhaps in Part I or Part II.
Question b) seems to be that you repeat the task like in a), for each X0=A+B*I, but a double-loop doing that was already illustrated here. You might try studying it some more.
acer

So, you want to save all the final (10th) iterates, for each A and B?
If that's true, then simply create an Array, outside the double loops. (Is it 101x101 in size?) Then after each k-loop finishes, save the final Xnew to the A,B coordinate of the Array. Then return that Array at the end of the procedure.
This Array would then contain the 10th and final Newton iterates, for all the X0=A+B*I initial points.
You could then do some sort of complex plot. Perhaps the idea was to show you how the various subregions of the A+B*I complex domain gets mapped to final iterates.
ps. For anyone who was shocked that the derivative of f is computed each time through the triple loop, sorry! Of course, it would be so much more efficient to get diff(f,x) just once, outside all the loops. For shame!
acer

I'm guessing that the idea is to start Newton's method using each of the points X0=A+B*I in the complex plane. So, the chances are greater to find a point that converges to a root. And perhaps several roots may be found.
N := proc(f)
local Xnew, X0, incX, A, B, k, results;
results := {};
# These two outermost loops are used to create the X0 initial points.
for A from 0 by 0.02 to 2 do
for B from -1 by 0.02 to 1 do
X0 := evalf(A + B*I);
# Now do ten Newton iterations, using X0
for k to 10 do
incX := evalf(eval(f, x = X0))/evalf(eval(diff(f, x), x = X0));
Xnew := X0 - incX;
X0 := Xnew;
end do;
# If it's "good enough", and not found already,
# then put the final iterate in the set of results.
if abs(incX) < 0.1000000000*10^(-6) and not member(evalf[7](Xnew), results) then
results := results union {evalf[7](Xnew)};
# Print which initial point converged to which new result.
printf("initial point %Ze converged to %Ze\n", evalf(A + B*I), Xnew);
end if;
end do;
end do;
# Return the results
return results;
end proc:
# try it out
N(sin(x));
Look, this may not be what you were assigned to do. It was unclear. If it's close, then study it, then when you understand it you should be able to alter it appropriately.
acer

The difference between the exact and floating-point eigenvector results, for say your 5x5 example, is not a condition number issue.
The reason for the differences you see in the exact vs floating-point results is simply that, for any given eigenvalue, the eigenvectors live in a linear subspace. Multiples of single eigenvectors, or linear combinations of eigenvectors in that relevant subspace (eigenspace) also happen to be valid eigenvectors for that given eigenvalue.
Take your 5x5 exact example. The exact eigenvalues are {3/4,0,23/20,1,1}. The eigenvalues {3/4,0,23/20} each have a single eigenvector associated with them, in the corresponding column of result Matrix e. You can scale each of those, and the result will still be an eigenvector of the same eigenvalue. The proof is easy, and falls out of the linear quality of the definition,
M * x = lambda * x
Similarly for the double eigenvalue {1,1}. It has two linearly independent eigenvectors associated with it (not always true, but true in this example). Any linear combination of such basis eigenvectors will produce another eigenvector in the same eigenspace associated with that same eigenvalue (ie, 1). This too falls right out of the linear quality of the definition above.
I don't think that numerical conditioning has much, if anything, to do with your 5x5 example when computed in floating-point.
I use Linux, and the particular scaling or linear combinations done on Windows might differ. But here is a sequence of elementary column operations with which I was able to transform the eigenvectors of evalf(L) into the floating-point equivalent of the results of the exact Matrix L.
Lf := evalf(L):
vf, ef := Eigenvectors(Lf):
Digits:=15: QQ:=ColumnOperation(simplify(ef),1,1/0.4472135955):
ColumnOperation(QQ,3,1/0.5920027785,inplace):
ColumnOperation(QQ,4,-1/(.7071067812),inplace):
ColumnOperation(QQ,2,-1/0.8164965809,inplace):
ColumnOperation(QQ,[5,2],0.03887187786,inplace):
ColumnOperation(QQ,5,(0.5/0.706304985706392),inplace):
ColumnOperation(QQ,[2,5],-1,inplace):
ColumnOperation(QQ,5,1/0.5,inplace):
ColumnOperation(QQ,[5,2],1,inplace):Digits:=10:map(fnormal,QQ);
Note that I only had to scale the 1st, 3rd, and 4th columns, relating to the eigenvalues of multiplicity one, to get an accord. To handle the eigenvalue 1, which has multiplicity two, I took linear combinations of just those the 2nd and 5th columns.
You may also enjoy playing with crude forward error estimates, using the definition, ie,
Norm( Lf . ef - ef . DiagonalMatrix(vf) );
You may also enjoy playing with the EigenConditionNumbers routine. These might help, if you doubt the floating-point calculations for your 50x50 case.
acer