If I understand what you want, here's one way to get it:

 

Use the following two procedures. The first is a help procedure for Q. The second plots the orbits by A of each the "points" L[1],L[2],...,L[k] in the list out to N times. That is it plots the points
L[1],A(L[1]),A^2(L[1],...,A^N(L[1]),
L[2],A(L[2]),A^2(L[2],...,A^N(L[2]),

#This is not as fast as Carl's but should work nicely for small matrices.

 

Orbit:=proc(A::Matrix,L::list,N::posint)
local P,j,Orb;
P:=Vector(L);
Orb:=L;
for j from 1 to N do
Orb:=Orb,convert(A^j.P,list)
od;
Orb;
end proc:

Q:=proc(A,L,N)
local O,i;
O:=NULL;
for i from 1 to nops(L) do
O:=O,Orbit(A,L[i],N);
od;
plots:-pointplot([O]);
end proc:

#Example:

 

A:=Matrix([[cos(Pi/8),-sin(Pi/8)],[sin(Pi/8),cos(Pi/8)]]);
L:=[[1,0],[2,0],[3,0]];

N:=10;

#Now see what happens when you execute the command

Q(A,L,N);

 

---Edwin

This procedure, I call Orbits,  has as input an nxn matrix A, a list v of length n, and a number N.

It returns the list [v,Av,A^2v, A^3v,...,A^Nv];

The problem is that starting with a list v to obtain the effect of A on v in list form one must write

convert(A.Vector(v),list).  Using lists rather than vectors makes it easier to plot the points.

Anyhow here's my procedure. (With more work one can replace the single list L by a list

of several points and get the orbits of all points in the list. Here L just represents a single point

in n-space if A is an nxn matrix. N + 1 is the number of points in the orbit. This counts the point L.

Orbits:=proc(A::Matrix,L::list,N::posint)
local P,j,Orb;
P:=Vector(L);
Orb:=L;
for j from 1 to N do
Orb:=Orb,convert(A^j.P,list)
od;
[Orb];
end proc:

 

Here is a very simple example

A:=Matrix([[2,0],[0,2]]):
L:=[1,1]:
Orbits(A,[1,1],5);

[[1, 1], [2, 2], [4, 4], [8, 8], [16, 16], [32, 32]]

I notice that if I set

Digits:=100;

then all methods give the same result to the 49th decimal place.

However for numerical integration I believe the correct command is

evalf(Int(f, x=a..b, ...))

This form gives the same answer for all three cases.

Edwin

I'm not quite sure what you want, but if I follow you correctly this might work:

 

with(LinearAlgebra):
n,m:=2,4;
M:=Matrix(n,m,(i,j)->2*i+j^2);  #just an example
A:=convert(M,listlist);
L:=map(s->add(j,j=s)/m,A);
V:=Vector(L);  #convert the list to a vector

;

Thanks for the suggestions, but I really had something more general in mind like the following (based on the use of evalindets which Joe Riel brought to my attention in a different context. I also use an idea similar to that suggested by John Fredsted Here's what I have now: Note that it works with sets of sets of lists of sets and all combinations of such: restart: Smp:=proc(A) evalindets(A, algebraic, x->trunc(x*10^10)); end proc: AreEqual:=proc(X,Y) evalb(Smp(X)=Smp(Y)); end proc: IsMemberOf:=proc(x,S) member(Smp(x),Smp(S)); end proc: Digits:=50: x:=evalf(Pi): x2:=x+Float(6.0,-49): AreEqual(x,x2); true S1:={{[x,evalf(sqrt(2))],[x2,evalf(sqrt(2))]}} union {{x} union {[1.0,0]}}: S2:={{[x2,evalf(sqrt(2))],[x2,evalf(sqrt(2))+Float(7.0,-30)]}} union {{x} union {[1.0+Float(7.0,-29),0+Float(1.0,-30)]}}: AreEqual(S1,S2); true Note evalf[10] doesn't work. The problem seems to lie with 0. I just raised this question in another thread. evalb(evalf[10](S1)=evalf[10](S2)); false