I'd select the text, copy it and then paste in into a text file such as code.txt. Then read it into maple using the read command, something like read(cat(currentdir(),"/code.txt")); (assumes your maple worksheet and code.txt are in the same directory). You may have to edit the text file slightly, e.g., what look like backquotes don't seem to be, so replace them with actual backquotes (according to the syntax error).

Also conformal in the plots package will do the grids in to new curves as you wanted. See also conformal3d.

After creating two plots, e.g. by

p1:=plot(x,x=0..1,color=blue);p2:=plot([[0.5,0.5]],style=point,symbol=solidcircle,symbolsize=30,color=red);

then plots:-display(p1,p2) plots p2 on top of p1, and plots:-display(p2,p1) plots p1 on top of p2.

Simplest way is

polarplot(sin(t)/cos(t)^2,t=-Pi/2..Pi/2, coordinateview=[0..2,-Pi..Pi]);

This may not be as efficient since it generates 256 matrices and then eliminates the zero matrices. But it does the conversion to the required notation at the same time.

 Full output of S is removed but starts
{E11, E12, E13, E14, E21, E22, E23, E24, E31, E32, E33, E34, E41, E42, E43, E44, E11+E12...

Download bits.mw

Edited: I realize I should have just gone from 1 to 15 (not 0 to 15) and then I woudn't have to remove the zero matrices.

The returned module pds has routines to plot the solution, e.g., pds:-plot(...). For details and examples see the pdsolve/numeric help page ?pdsolve,numeric

Most of Maple (including the plot3d command) uses the "Math" convention as on the ?coords help page, i.e., the azimuthal angle second.

But the VectorCalculus package uses the "Physics" convention as on the ?VectorCalculus,Coordinates help page, i.e., azimuthal angle third. I believe the Physics package also uses this convention, though a quick glance at the help pages did not clarify this.

This is a definite problem for Maple, and I have seen some commentary about improving the situation, though I don't recall where. A first step would obviously be to put some warnings on each of these help pages!

I used 

read "stanform";

and got some immediate errors because Maple's ditto operator used to be " but is now %. I didn't attempt to replace the " with %, but it might be possible to make it work if you are familiar with Maple. In general old Maple code still runs, but there are a lot of lines of code there and I didn't look in detail through it.

According to the help page ?pdsolve,numeric, the default method is used unless you specify one of the others (those in your list). Under "Details" on that page is some more information:

"The default method uses a second order (in space and time) centered, implicit finite difference scheme to obtain the computed solution. The number of points in the stencil of the finite difference scheme is one greater than the order of each equation."

Some references are give at the bottom of the page, though they seem to be rather general.

PDETools:-dcoeffs(lhs(eq),q(T));

gives

324.6463527, 3.589858529*10^12, 4.012505275*10^11

eval substitutes N=infinity in each term, and concludes each term and the sum is zero. Limit is smarter:

limit(G(N),N=infinity); returns without an answer. I tried some assumptions here withouth success. On the other hand, you can do it numerically and approximately:

evalf(eval(G(N),N=10000)) gets it right to 3 decimal places. 

For 8x8

InvM:=Adjoint(M)/Determinant(M):

returns quickly, but then InvM[1,1] returns with "[Length of output exceeds limit of 1000000]"

So I think the answer is too complicated to do anything with, no matter how much memory.

On the other hand,

simplify(InvM.M);

gives the identity in about 150 s on my laptop.

You need to make an actual second point using Point1:=disk( ...  , color =black), then rotate it and display it with Point2. Rotate doesn't get rid of the color of the object, and the display color won't override a specified color in a plot object.

David.

In Maple 16:

restart;Digits;
10
evalf(Sum(I^n*BesselJ(n,2*sqrt(34))*exp(-I*n*arctan(5/3))*exp(I*n*Pi/4),n=1..30));
evalf(Sum(I^n*BesselJ(n,2*sqrt(34))*exp(-I*n*arctan(5/3))*exp(I*n*Pi/4),n=1..100));
evalf(Sum(I^n*BesselJ(n,2*sqrt(34))*exp(-I*n*arctan(5/3))*exp(I*n*Pi/4),n=1..1000));
0.05400552045 - 0.7873756700 I
0.05400552043 - 0.7873756700 I
0.05400552042 - 0.7873756700 I

L:=[{a},{b,c},{d,a,f},{b,d}] ;

look:=L[1][]:   #thing to look for
for i from 2 to nops(L) do
 if has(L[i],look) then print(L[i] minus {look}); break; end if;
end do;

gives {d,f}

which I think is what you wanted.