Another workaround is to set Digits:= 20.

Assumptions via assume on one variable are not cumulative. Your second assume command erases the effect of the first. The second one is unnecessary anyway: It is implied by the first. If you actually want to make multiple assumptions on a single variable, put them together in one assume command or use the additionally command. See the examples at ?assume for details. If you remove the second assume, you will get your expected result.

By the way, you misspelled alpha as alph, but that doesn't change the result.

I agree with Preben about assuming. I think that it makes for cleaner, easier-to-read, less-error-prone, and more-flexible worksheets. Although there is still a place for assume if you're trying to get Maple to do a symbolic proof.

This is a bug in the Standard GUI. Your commands run fine in the command-line Maple. If you have access to the Classic GUI (I don't), try that also.

A warning to anyone trying to investigate this: If you try to probe this, say with printlevel:= 100, the GUI will lose connection with the kernel (without any message), and you'll be stuck. Setting interface(prettyprint= 0) does not help. Another curiosity: The output, i.e. the returned value, of the OP's commands is literally the symbol `Non-fatal error while reading data from kernel.`; that is what soltn is set to; it is not an error message in the usual sense. Therefore, it must the kernel that is sending the message!

There are two independent variables, t and x. So they need to be 3D plots. I don't think two 3D plots, together, of functions that are so close will be very useful. But it can be done if that's what you really want. You can have two side-by-side 3D plots on different axes. Another option is a plot of the difference between the two functions.

There are two immediate problems:

First, the initial condition P(0) = P uses P in two different ways. You need to use a different letter.

Second, the exponential function is represented in Maple as exp(...), not e^(...).  If you correct this to exp(-1.15*10^12), the resulting number is so incredibly small that it is difficult (not impossible) to work with. So, before we go further, I have to ask if this is what you really meant.

As a statement, to return from a procedure:

if `and`(0<z, abs(y)<=6, 0<x,x<1) then return else return 0 fi

As an expression, which may or may not appear after return:

`if`(`and`(0<z, abs(y)<=6, 0<x,x<1), NULL, 0)

NULL is equivalent to [][]; I prefer the latter, which is easier to type and doesn't clutter the code with words in ALL CAPS. In several contexts in Maple, a NULL is implied by the absence of anything else; such is the case with the first return statement above. The if-statement bracket fi is if backwards, and is equivalent to end if. And `and`(cond1, cond2, ..., condN) is equivalent to cond1 and cond2 and ... and condN, with the former being less to type when there are more than two conditions.

Another way --- one of many.

restart;
r:= [3,3,3,3,4,4,4,3,3,3,3,3,3,3,3,2,4,3,1]:
T:= ListTools:-Classify(nops, [ListTools:-Categorize](`=`, r));

T := TABLE([1 = {[1], [2]}, 4 = {[4, 4, 4, 4]}, 13 = {[3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3]}])

m:= max(indices(T));

m := 13

op ~ (T[m])[];

3

I wrote the above to cleanly handle the multi-modal case (as Kitonum also did), which adds a bit of complexity.

Yes, dsolve with type= series says that the positive-term series that you posted is the solution. I think that I found the sign bug in your code.  You have a line

u[m,k+1]:= -R[m](seq(u[i,k],i=1..nEq))-A[m,k];

I think that the last term should be +A[m,k]. If I make this change, then I get the positive-term series that you want.

Maple is saying that there's a singularity near x = 2.5. Do you have reason to doubt that? Here's the code and plot. You can see that it's approaching an asymptote.

restart;
dq:= diff(y(x),x) = x^2 + y(x)^2 - 1:
DEtools[DEplot](
    dq, y(x), x= 0..2.5, {y(0)=0}
   ,view= [default, -1..5]
   ,arrows= line
   ,title= 'Slope*Field'
);

Just apply coeff twice:

S:= A*sin(x)*cos(z) + B*sin(2*x)*cos(z) + C*sin(3*x)*cos(2*z):

coeff(coeff(S, sin(2*x)), cos(z));

                                                                 B

The command that you want is ?subsindets. For example,

subsindets(expr, anything, simplify);

As you're learning to use this, you will quickly come to a situation where you don't want to apply the transformer to absolutely every subexpression, yet there is no convenient type (the type in the example is anything) that restricts it to what you do want. At that point you'll need to learn ?frontend.

Both subsindets and frontend are built-in kernel commands, so they are more efficient than user-defined procedures.

If you literally want to apply map to all levels of an expression, you can easily make a recursive version of map. For example,

DeepMap:= proc(f, expr)
    option system, remember;
    f(`if`(expr::atomic, expr, map2(procname, args)), args[3..-1])
end proc;

(This one applies f to the outermost level also.) (This handles the basic map, not any of map's indexed versions.)

The subs command only does "syntactic" substitutions, i.e. the item being substitued for must occur as a distinct entity within Maple's internal representation of the expression. It is more of a "programmer level" command than "user level", and it makes no claim as to the mathematical validity of what it does. To perform algebraically valid substitutions, use algsubs.

Use ?LinearAlgebra,LeastSquares:

   X:= LinearAlgebra:-LeastSquares(A,Y);

It is a nuisance that geom3d requires you to label every object. The object model used in that package is much older than the modern one based on modules. In the modern module model, an object "is what it is", so to speak, regardless of what, or even if, you name it; in the older model, an object's properties are based on looking it up by name. (To "deconstruct" it, it seems that the evolution of the object model is analogous to the evolution of social classes.) But you can reuse the labels in geom3d, so the easiest way to do what you want is

    [seq](sort(Equation(TangentPlane(P, S, point(A, pt[])), [x,y,z])), pt in L);

(where all points are named A and all planes are named P).

If you do require distinct names for each object (so that you can look up their properties without recreating them), it is possible to construct names based on the numerical values of the point using nprintf:

[seq](
       sort(
            Equation(TangentPlane(nprintf("P%a",pt), S, point(nprintf("A%a", pt), pt[])), [x,y,z])
       )
      ,pt in L
 );

Unique names can also be created with cat, which is what John May did.

I made two very small changes to your procedure, and I almost got the series that you indicated, the only difference being that I got an alternating series. Here are a few tips to help you debug your procedures:

1. Indent after each level of do loop.
2. Use the trace command -- in this case trace(Ad) -- before running the procedure. This will show you the result of each assignment (":=") statement. 

And here are the small changes that I made:

1. Do not use unevaluation quotes; for example use u[i,k] rather than 'u[i,k]'.
2. Use add, not sum.