@Carl Love I like your revisions quite a bit.  Thanks!

@Kitonum Is there a reason to believe that the m vectors returned by Basis() contains the original n rows of the matrix?

@sideshow Your code attempts to solve the Poisson problem by mapping a semicircular domain into a rectangle through the usual Polar <-> Cartesian change of variables. You should know, however, that after such a mapping, the Laplacian operator is no longer u_rr + u_θθ. It changes to u_rr + (1/r) u_r + (1/r^2) u_θθ; see, e.g., Wikipedia.

If you want to push in that direction, you will have to change your finite difference scheme to handle this modified differential operator but I wouldn't advise it. The finite difference method works well for simple PDEs in simple domains. The finite element method was invented to handle more complex PDEs and geometries.

@sideshow The purpose of algo_12_1() is to solve a Poisson problem on a recangle [a,b]x[c,d], and the boundary condition g is expected to be the values of the solution along that rectangle's edges.  You can't use it to solve a problem on a semicircular domain.

The error message you are seeing arises because the procedure attempts to evaluate your g on the rectangle's boundary, but g returns unevaluated there, by design.

@Carl Love You are correct in focusing on the 10-digit software-float arithmetic versus hardware-float arithmetic.  The difference is of significance particularly in the sample problem solved (from the book) which sets TOL = 1e-10 for the stopping criterion, therefore necessitating higher than 10-digit precision to get meaningful results.

But that's not the only source of difficulty in OP's code.  There we have several stanzas of the type

if abs(w[i,m-1] - z) > NORM then
       NORM := abs(w[i,m-1] - z);
       w[i,m-1] := z;
 end if;

which should have been

if abs(w[i,m-1] - z) > NORM then
       NORM := abs(w[i,m-1] - z);
 end if;
 w[i,m-1] := z;

There may be other errors in that code, but as I noted earlier, I have not really attempted to debug it.

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