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.

When applied to an rtable (Matrix, Vector, Array), nops is not the number of elements; rather, it is the number of operands in the internal structure. You can see the internal structure of A via op(A). For the number of elements, use rtable_num_elems(A, 'All').

I'd like to show you the relevant commands, but I also want you to do your homework. So I'll compromise by posting the following without comments, most output elided, and without an attached file.

restart;
X:= t-> 3*cos(t):  Y:= t-> 2*sin(t):  t0:= 0.5:
kappa:= eval(VectorCalculus:-Curvature(<X(t),Y(t)>), t= t0):
dist:= (X1,X2)-> sqrt((X1[1]-X2[1])^2+(X1[2]-X2[2])^2):
slope:= (D(Y)/D(X))(t0):
x0:= X(t0):  y0:= Y(t0):
eq1:= (Cy-y0)/(Cx-x0) = -1/slope:
eq2:= dist([Cx,Cy],[x0,y0]) = 1/kappa:
fsolve({eq1,eq2}, {Cy,Cx});
                        {Cx = 4.13904333495107, Cy = 2.19319067267216}
fsolve({eq1,eq2}, {Cx,Cy}, avoid= %);
                        {Cx = 1.1264520363911743, Cy = -0.27548851825535093}
assign(%);
tan_eq:= slope*(x-x0)+y0:
Ellipse:= plot([X,Y,0..2*Pi], color= blue, scaling= constrained, thickness= 2):
Circle:= plot([Cx+cos(t)/kappa, Cy+sin(t)/kappa, t= 0..2*Pi], color= red):
Center:= plot([[Cx,Cy]], style= point, symbolsize= 25, symbol= cross, color= red):
Tangent:= plot(tan_eq, x= 0..4, color= green, thickness= 2):
plots:-display([Ellipse, Circle, Center, Tangent]);

There are some problems with your terminology. You say gg is a vector. It doesn't make sense to raise a vector to a power. Perhaps g is a vector and gg represents the vector's dot product with itself, which would be a scalar. What you call the "gradient of psi" is the derivative of V with respect to psi. You can get this as

diff(V, psi);

But that will not express the derivative in terms of V.

You expression GVpsi is off by a factor of V. In other words, the numerator of the coefficient in front of the expression should be V instead of 1. It is difficult to get Maple to produce this expression as the derivative of V. However, you can verify it by subtracting diff(V, psi) and simplifying, obtaining 0. The expression can be derived "by hand" by a process called "logarithmic differentiation", which you can look up in a calculus textbook, which is often used to take the derivative of complicated expressions of the form f(x)^g(x).

Download Psi.mw

All algebraic grouping in Maple, no matter how deeply nested, is done with ordinary round parentheses (), not square brackets [] or curly braces {}.

Trying to separate assigned and unassigned names can get very confusing, especially if you try to generalize it by putting the names in a set and then iterating over that. Here's a much simpler method to set certain SpecialNames to 0 in an expression ex if they are not assigned.

SpecialNames:= {h||(251..256)}:
evalindets(ex, name, x-> `if`(x in SpecialNames, 0, x));

You can make up weird situations where this won't work, like if one name in the set is assigned to another name in the set. But you're not going to do that weird thing, right?

Let me know how that works for you.

chrem([3, 3, 10, 0], [7, 11, 13, 7+11+13]);
                              465

P1:= plots:-polarplot([1,t,t=3*Pi/16..7*Pi/16]):
plots:-display([
     P1,
     plots:-polygonplot(
          < < 0 | 0 >, plottools:-getdata(P1)[-1], < 0 | 0 > >,
          colour= violet
     )
]);

See command ?assigned.

See ?LinearAlgebra,QRDecomposition.

If you like the way code is entered in Classic Worksheet (I do), then you can select the style "Maple Input" from Standard Worksheet. I also like Code Edit Regions (select from the Insert menu) for sections of code between 5 and 100 lines.

Expanding on my earlier reply, you could take all the usused elements of set X and assign them to NULL or some other default value like the empty set {}.

T:= table([op(2, eval(T))[], (X minus {indices}(T,'nolist') =~ ''{}'')[]]);

Note that that's two nested pairs of single quotes, not one pair of double quotes.