I think that what you want is

plots:-spacecurve([seq]([x[k],y[k],z[k]], k= 1..200));

But I don't know if that's "in the three-dimensional phase domain", and it's certainly not a "phase plane".

Anyway, please let me know if the spacecurve is what you want.

You need to make some assumptions for the IsDefinite to work.

assume(w1 >= 0, w2 >= 0, w3 >= 0, w4 >= 0, w5 >= 0, w6 >= 0);

I don't know if the overall program will work after these assumptions are made (you didn't attach the uploaded file), but that's the first hurdle.

Note that L(w1,w2,w3,w4,w5,w6,lam) is evaluated before being passed to NLPSolve. Thus, the fact that you did specify the nonnegativity assumptions in the call to NLPSolve does not cause those assumptions to be used while evaluating the arguments to NLPSolve.

This is a bit tricky; perhaps not suitable for the Maple beginner.

F:= (N::evaln(list(`=`)))-> (convert(N, `local`)@lhs = rhs) ~ (eval(N)):

u:= [a=x, b=w]:
p:= [1=4, 3=2, 5=2]:
F(u), F(p);
      [u(a) = 1, u(b) = w], [p(1) = 4, p(3) = 2, p(5) = 2]

seq(Matrix([[seq(a[k], k = x+m .. y+m)]]), m = 0 .. 2);

^^^^^^^

Obviously, the result of such a command is a sequence of Matrices, not a Matrix.

Here's my shortest code:

< seq(< a[x+m..y+m] >^%T, m= 0..2 ) >;

In particular, note the for any list L and valid indices a and b,

[seq(L[k], k= a..b)]

is the same thing as simply

L[a..b]

In Maple, when a negative real number is raised to a fractional power with an odd denominator, the result is not real. As it says at ?arithop or ?^,

For non-integer (complex) numeric powers y and (complex) numeric x, the expression x^y is evaluated as exp(y*ln(x)), where ln(x) is evaluated using the principal branch of the logarithm.

The function surd can be used to compute real odd roots of negative real numbers.  See ?surd for more information.

Since ln(-1) = I*Pi, the upshot of that first sentence is that for integer n, (-1)^(1/n) = cos(Pi/n) + I*sin(Pi/n):

radsimp((-1)^(1/3));
                         1   1    (1/2)
                         - + - I 3     
                         2   2      

cos(Pi/3)+I*sin(Pi/3);
                         1   1    (1/2)
                         - + - I 3     
                         2   2         
But

surd(-1, 3);
                               -1

or

use RealDomain in (-1)^(1/3) end use;
                               -1

(see ?RealDomain).

In this case, you need to circumvent the automatic simplification of x^0 to 1. You do this with

coeff(b, x, 0).

We make Points1 into an (M+1) x (N+1) x 3 listlistlist:

A:= [seq](Points1[(M+1)*k+1..(M+1)*(k+1)], k= 0..N);

plots:-surfdata(A);

And you can add the same options that you used for the pointplot3d command.

Please remove the with(linalg) from your worksheet. It works without it.

If they can be considered the same thing and also be considered the same as 1, then it is trivial:

Re ~ (sol);

Then proceed with whatever else you were doing with the repetitions.

If your goal is simply to remove the repetitions from a list L, and you don't care about the order, then it is trivial: convert to a set:

{L[]};

If you want to keep it a list, but you still don't care about the order, it is still trivial: convert to a set, then convert back to a list:

[{L[]}[]];

When order doesn't matter, removing repetitions, whether or not they exist, is much easier than detecting or finding them: The former can be done entirely in Maple's kernel.

An rtable is a Matrix, Vector, or Array. A list is not an rtable. FindRepetitions returns a list (so converting it to a list doesn't do anything). Why don't you use the numelems command that acer recommended? If your Maple version is less than 15, use nops for a list.

If the implicit-form solution can be solved for the independent variable, as it can in this case, then there are several more very easy ways to get it.  Here's three: (Sorry, explicit uploading of worksheets seems broken right now, so I am uploading this in plaintext (prettyprint=1).)

restart;
ode:= x - y(x)/diff(y(x), x) = y(x);
                            y(x)         
                      x - -------- = y(x)
                           d             
                          --- y(x)       
                           dx            
(1)
dsolve(ode);
                       /        /     x    \      \
             y(x) = exp|LambertW|- --------| + _C1|
                       \        \  exp(_C1)/      /

isolate(subs(y(x)= y, %), x);
                     x = -y (-_C1 + ln(y))

(2) The name of the conversion below is y_x regardless of the names actually used in the ODE.

convert(ode, y_x);
                       d              x(y)
                      --- x(y) = -1 + ----
                       dy              y  
dsolve(%);
                    x(y) = (-ln(y) + _C1) y

(3) The order of  the substitutions below is significant.

subs(diff(y(x),x)= 1/diff(X(y),y), y(x)= y, x= x(y), X= x,  ode);
                             / d      \    
                    x(y) - y |--- x(y)| = y
                             \ dy     /    
dsolve(%);
                    x(y) = (-ln(y) + _C1) y

Download implicit_dsolve.mw

For any complex y, y*conjugate(y) = |y|^2 (although I would've guessed that you already knew this identity). The identity can be applied automatically via

simplify(z);

The resulting equation can be solved for abs(y) via

solve(%, abs(y));

indicating that there is a circle of solutions in the complex plane.

How about this? With the data in your worksheet, this took exactly 200 iterations. The resulting plot is a straight line (as well as can be seen with the naked eye). Note the the default for Norm is the infinity norm, which is equivalent to checking the absolute values of all of the entries. In other words, the convergence criterion is that the distance for both coordinates of all n-1 points is less than epsilon.

ErrX:= Vector(n, [0,1]):
ErrY:= Vector(n, [0,1]):
rho:= 1.2:  #relaxation factor
epsilon:= 10.^(-5):
while LinearAlgebra:-Norm(< ErrX, ErrY >) >= epsilon do
    for k from 2 to n do
        Xk:= X[k] - rho*GX(X, Y, k)/GXX(X, Y, k);
        ErrX[k]:= Xk - X[k];
        X[k] := Xk;
        Yk:= Y[k] - rho*GY(X, Y, k)/GYY(X, Y, k);
        ErrY[k]:= Yk - Y[k];
        Y[k]:= Yk
    end do
end do;

restart;
m:= proc(r::nonnegint, c::nonnegint)
   option remember;
   if r*c=0 then 0 else c*thisproc(r-1,c)+a*thisproc(r-1,c-1) end if
end proc:
m(0,0):= 1:

M:= Matrix(5,5, m);

Matrix indices start at 1, not 0, so the (trivial) 0th row and column are not shown. If you want indices that start somewhere other thn 1, you need to use an Array. But then you won't get the neat tabular layout in the GUI.

Your procedure myp can be reduced to

myp:= (a,b,c)-> zip(`+`, a, b)*c;

Is it possible to use something like
#map2(myp,[$1..4],{[x,y,z,w],[x2,y2,z2,w2],[x3,y3,z3,w3]},[4,10,0.5]);
to achieve this?

seq(map2(myp, [$1..4], {[x,y,z,w],[x2,y2,z2,w2],[x3,y3,z3,w3]}, s), s= [4,10,0.5]);

(I don't know why the line above breaks bad. It appears as one line in the editor, and I've tried several times to fix it.)

Is there a document I can read apart from the ?map to look for some general tips to 'simutanous' computing both symbolically and numerically in Maple?

See ?operators,elementwise and ?curry.

And is there any documentation on using multicore CPUs?

See ?Threads.