@Carl Love Thanks for the reference.  It had not occurred to me to make a connection between this and pseudoinverses.  I have learned somthing new today.

And as @vv has pointed out, although the Wikipedia page does the Ax=0 case, extending it to the Ax=b  is straightforward.

@Bendesarts Try the arrows=none option with DEplot().  See the "Starbucks Girl" for an example.

@Bendesarts That message says:

Error, (in dsolve/numeric) x(t) and x cannot both appear in the given ODE

This is an error in your code that Carl pointed out in your initial post, then I pointed out, but it has not been corrected yet.  You should be more careful with the details. Sometimes you write x, and sometimes you write x(t).  That's not good.  They should all be x(t).

@Bendesarts The RealDomain command is not useful in this case.  I used it in the previous case to locate the equilibrium points which were five of them.  But this other equation has only the origin as the equilibrium, so there is nothing to calculate.

As to initial conditions, you want to draw multple orbits that start at the plotting rectangle's boundary.  To generate points along the left and right edges we do:

ic := seq([x(0)=xmax, z(0)=h], h=-zmax..zmax, 0.1),
        seq([x(0)=-xmax, z(0)=h], h=-zmax..zmax, 0.1);

I will let you add points along the top and bottom edges.  This will fill the region outside the ellipse.  See what you can do to fill the inside, and if you need help, come back and ask.

As to adding points on an oribit, that is beyond what DEplot can do.  To add a point you need to actually solve the differential equation, like this:

dsol := dsolve([EqSys[], x(0)=0.6, z(0)=zmax], [x(t), z(t)], numeric);

Let's pick three points along the resulting orbit corresponding to t=0, t=1, and t=2.  We call them p0, p1, p2:

p0 := eval([x(t),z(t)], dsol(0));
p1 := eval([x(t),z(t)], dsol(1));
p2 := eval([x(t),z(t)], dsol(2));

Plot these points by calling the pointplot() function.  That function lives in the plots package, so we load it first:

with(plots):
pointplot([p0, p1, p2],
    color=[red,green,blue], symbol=solidcircle, symbolsize=20);   plt2 := %:

The plt2 := %: at the end of that line names the plot of points plt2.   The plot produced DEplot we call plt1:

DEplot( ... as before ...);  plt1 := %:

Now put plt1 and plt2 togheter:

display([plt1, plt2]);

Done!

@Bendesarts This most recent version of your equations still contains many of the syntax errors that were pointed out by Carl.  Pay attention!  Details matter!

After you have fixed those errors, you will find out that the graph produced by Maple is nowhere near what you have shown.  The reason is that there is a sign error in your eqz.  The plus sign in front of w should be a minus.

Here is the fixed version:

restart;
r := sqrt((x(t)/a)^2+(z(t)/b)^2);
eqx := diff(x(t),t)=alpha*(1-r^2)*x(t)+w*a/b*z(t);
eqz := you do it
params := alpha=1, beta=1, a=0.4, b=0.2, w=1;
EqSys := eval([eqx,eqz], [params]);
xmax := 0.5;  zmax := zmax; tmax := 8;
ic := [x(0)=0.5, z(0)=zmax], [x(0)=0.01, z(0)=0];           # see the note below
DEtools[DEplot](EqSys, [x(t),z(t)], t= 0..tmax, [ic],
    linecolor=black, thickness=1, stepsize=tmax/(3*48),
    x(t)=-xmax..xmax, z(t)=-zmax..zmax, scaling=constrained);

Note: To specify the initial conditions systematically, see my "Starbucks girl" post.

@Markiyan Hirnyk The value of x0 is immaterial.  In fact, just let x0 be 1:

z := x^2/(1 + (1-x^2)^2)^(3/2);

int(z, z=0..infinity);

simplify(%);

∞

This is in Maple 2015.1.

@nhungnhung 

restart;
with(Student[Calculus1]):
SurfaceOfRevolution(x^2+x-1, x=-3..1, axis=vertical,
    output=plot, transparency=0.4, axes=boxed);

"It does not work" does not tell me what it is that you see.  Can you be more specific about what goes wrong?

It would help if you explained what it is that you want to achieve with those calculations because as written, much of it makes no sense.

1. You define a function H and claim that it is constant.  How do you know that?  In fact, it isn't because omega is not constant.
2. You plot what appear to be orbits in the the q-p phase space, but they are not because omega is not constant.
3. You are calcuating the area under the supposed orbits.  Why do you need the area?
4. The meaning and significance of t* is unclear.
5. You should not assign to the letter I in Maple because I in Maple stands for sqrt(-1).

The main source of confusion is in treating a non-autonomous equation as if it were autonomous.  Some of the issues noted above arise due to that confusion.  The techniques you are attempting are not applicable to non-autonomous equations.

If you are not familiar with the terms autonomous and non-autonomous, you may safely think of an autonomous equation as one whose coefficients do not depend on time.  Your differential equation is non-autonomous since the coefficient omega depends on time.

@snowww Your "(now corrected)" version makes mathematical sense, however I don't trust it since it is dimensionally inconsistent.

The differential equation implies that F and K have dimensions of 1/x.  But if so, then the condition

    1 = C - F*dC/dx     at x=0

is dimensionally inconsistent.  I suspect that what you want is the dimensionally correct

    dC/dx = F*(C-1)     at x=0.

@snowww Boundary conditions specify conditions at boundaries.  What is x doing in DC(0)=F(C(x)-1)?

The expression F(C-1) in the boundary condition makes no sense.  Probably you meant something else.

@der Uwe I haven't used the VectorCalculus package but I would like learn about it.  As a practice problem I attempted to express the Navier-Stokes equations in cylindrical (not torroidal) coordinates in Maple but was unable to derive the well-known answer.  I would be interested in seeing how you did your calculations so that I may learn something from it.  Would you mind posting your worksheet?

@Robert Israel Thanks for confirming and reporting this bug.