@Christopher2222 As the text in the picture that you have posted says, that contour plot is in a rotating coordinate system.  The equations that you are plotting don't account for the rotation at all.  Without rotation there are no inertial forces, and without inertial forces there are no Lagrangian points.

NomNomPancake Your equations account for the combined gravitational potentials of the two point masses only.  Those will be adequate if the masses were held in place without moving, say by god.

In reality you know that if he (she?) lets go of the masses, they will begin to move.  If the initial conditions are right, then jupiter will begin orbiting the sun.  If a particle is placed in a Lagrangian point of the sun and jupiter system and rotates with them, it will feel not only their gravitational fields, but it will also be subjected to inertial forces (centrifugal and Coriolis).  Those are missing in your equations.  There won't be Lagrangian points without accounting for the inertial forces.

@Johnluo I assume that the m:x->... is a typo and it should be m:=x->...

After fixing that, change the last command to:

intsolve(myeq,u(x));

Within a few second, Maple returns the obvious solution u(x) = 0.  But you don't need Maple to tall you that.  This was pointed out in an earlier message bytomleslie.

@Josolumoh To make the assumptions available throughout a worksheet, we do:

restart;
assume(a > -1, lambda < 0, beta > 0);
int(t^a*exp(-beta*t)/(-lambda*t+1), t = 0 .. infinity);

and we get:

The notation a~ indicates that an assumption has been made on the symbol a.

@Josolumoh One does not "expand" something like "GAMMA(-a, -beta/lambda)" in the same way that one does not "expand"  sin(a/beta).  One treats it as a known quantity.  To get the numeric value of  GAMMA(-3, 5) we do evalf(GAMMA(-3, 5)) and we get 0.0000062638. 

If you want to explore more, the precise definition of GAMMA(a, z) is given in Maple's help page on GAMMA.  You will see that it is defined in terms of the hypergeometric functions.  In special cases the latter are expressible in terms of the so called Exoinential Integral functions, that is the Ei(4,5) that you have noted.  Look up Maple's help on Ei to see what it is.  Normally you wouldn't want to be bothered with these details -- to get numerical values you just apply evalf as I have noted above.

@Maple4evergr8 This may depend on your operating system.  I use Ubuntu Linux in which plotsetup(window) is accepted but when attempting to plot, it fails with "no plot device driver for plotdevice=window".  To list all possible plot devices, we do:

plotsetup(help);

In Ubuntu I get a large list of devices, among which are:

    xwindow, x11, X11

which so far as I can tell are synoyms.  Therefore in a terminal window we do:

plotsetip(x11);
plots:-animate(plot, [a*x^2, x=-1..1], a=-1..1, frames=40);

and get a window pop up with a nice animation.

@Amal You haven't defined C[13].  What is it? 

You wrote: "I do not get the right result".

Why do you think that is not the right result?

@acer That's very clever.  Here is a thumbs up!

On computing the equilibria you get:

The first of these says that  (0,y,0) for any y is an equilibrium, that is, you have a continuum of equilibria.  In the rest of the worksheet you proceed as if there are two equilibria only, and treat the "y" in pts as if it were a number, which it is not.  Near the end you have:

indicating again that you think you have two equilibria.  You need to clarify this for yourself.

@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!