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.

Construct a graph of 12 vertices, one for each set. I'll label the vertices simply 1 - 12 rather than S[1] - S[12]. Two vertices will share an edge iff the corresponding sets do not intersect. Then the answer is simply the size of the largest clique in the graph.

LL:=[

   [2, 4, 5, 8, 10], [1, 3, 6, 7, 11], [2, 5, 6, 9, 10], [1, 5, 6, 9, 11], [1, 3, 4, 7, 12], [2, 3, 4, 8, 12]

 , [2, 5, 10, 11, 12], [1, 6, 10, 11, 12] , [3, 4, 10, 11, 12], [1, 3, 7, 8, 9], [2, 4, 7, 8, 9], [5, 6, 7, 8, 9]

]:

use GraphTheory in CliqueNumber(Graph({seq}(seq({i,j}, j in LL[i]), i= 1..12))) end use;

                                                                 3

Looking at the graph with DrawGraph, one can easily see many 3-cliques. (3-cliques are easy to see because they are just triangles.) One example is {1,4,5}.

One problem in your code is that your comments block some of the essential code. In any line containing a #, anything to the right of the # is ignored. (So there's also no need to use quotation marks in your comments.) For example, you have a line

#` Step 4`  for i from 2 to M2 do

So, the part for i from 2 to M2 do will be ignored.

I almost never use 2D-input. I only use it when it is part of someone else's worksheet that I'm playing around with. The only way that I know how to get a new line is Control-Enter, same as in 1D input. When I use that, I don't get the extra space. What is the usual way to get a new line, the way that does give the extra space?

I've modified the above so that the axes don't change due to the erasure, and the original axes will be restored if the VIEW is reverted to default by a subsequent display. So it is now more like erasing part of a blackboard (and circles will remain circles), whereas before it was like cutting off part of the blackboard. I do this by keeping track of the extreme upper, lower, left, and right points as they pass through the transformer procedure. Then I add an invisible (color= white) plot of the two extreme points: lower left and upper right.

It probably goes without saying that this procedure should work for any 2D plot, not just BubblePlot.

Eraser:= proc(Plot::specfunc(function,PLOT))
# Uses the VIEW specification of a PLOT structure to erase points outside the view.
# Input: PLOT structure with VIEW; Output: The PLOT structure, modified.
   local
      x1,x2,y1,y2,newPlot
     ,minx:= infinity, miny:= infinity, maxx:= -infinity, maxy:= -infinity
     ,View:= indets(Plot, specfunc(anything, VIEW))
   ;
   if View = {} then  return Plot  end if;
   if nops(View) > 1 then  error "Multiple VIEWs found: %1", View  end if;
   View:= View[];
   newPlot:=
      plottools:-transform(
         subs(
            [x1,x2,y1,y2] =~ [op(1..4, map(op, subs(DEFAULT= -infinity..infinity, View)))]
           ,proc(x,y)
               if not(x::numeric and y::numeric) then  return [args][1..2]  end if;
               if x < minx then  minx:= x  end if;  if x > maxx then  maxx:= x  end if;
               if y < miny then  miny:= y  end if;  if y > maxy then  maxy:= y  end if;
               `if`(`or`(x1 > x,x > x2, y1 > y,y > y2), [undefined$2], [x,y])                  
            end proc
         )
      )(Plot)
   ;
   # Add two invisible (white) points at extreme corners of original plot so axes don't change.
   # plottools:-transform removes VIEW, so put it back.
   plots:-display(
       [plots:-pointplot([[minx,miny],[maxx,maxy]], color= white), newPlot]
      ,view= [op(1..2, View)]
    )
end proc: