It seems like you don't expend much effort trying to find errors on your own, because you have so many, and some of them are so obvious (such as unbalanced quotation marks, or exchanging P and p (refering to your previous Question)). You just complain here that your code "don't work". Well, it's not like we here at MaplePrimes have some automatic way to find the errors. We just apply the logical rules by hand, as you could learn to do.

I think that this accomplishes what you want:

plots:-display(
   seq(
      seq(
         [
            plottools:-polygon(
               [[i,j],[i,j+1],[i+1,j+1],[i+1,j]], 
               color= `if`(j::odd, magenta, white)
            ), 
            plots:-textplot([i+.5, j+.5, sprintf("%d",i*j)])
         ][],
         i= 1..10
      ), 
      j= 1..10
   ),
   axes= none
);

If instead you want the numbers to increase as you go down, then change the j range from 1..10 to -1..-10, -1 and change i*j to abs(i*j) (or -i*j).

Your code executes without error for me. The file is plaintext, and is easily readable in most browsers or editors; so I can assure you that it contains exactly one matrix without any syntax issues. The first three positive roots of the determinant are

[0.269373966731148794754757341200462674642148338669373608296012774382885541088668300033662376394882171035784896969388875153559058382051590463209287590804e-1, 0.990674477020063147446983680621401355331598629110322989492668885141230642823201150297702439572643428277168426737360736443335299618942000792413816520383e-1, .101120616665296170767949963559947228275460392424861673345779263832172742662163268222684076145562969014514409071557800322470075818964100491423123404014]

These results agree with VV's above. His technique is more robust in terms of round-off error.

Obviously the number of cycles of length k is closely related to the number of walks of length k that end at their own starting point. Clearly, this number of walks is the trace of the k^th power of the adjacency matrix. It only remains to remove the walks that aren't cycles and to divide out the 2*k permutations of each cycle. This formula does that:

#Number of k-cycles for a given adjacency matrix A
kCycles:= (A::Matrix, n::posint, k::And(posint, Not({1,2})))->
   add(
      (-1)^(k-i)*binomial(n-i,n-k) * 
         add(LinearAlgebra:-Trace(A[S,S]^k), S= combinat:-choose(n,i)), 
      i= 2..k
   )/2/k
:

Continuing with the example of the octahedron graph:

A:= GraphTheory:-AdjacencyMatrix(GraphTheory:-SpecialGraphs:-OctahedronGraph()):
n:= 6: #number of vertices 
seq(kCycles(A, n, k), k= 3..6);
                         8, 15, 24, 16
add([%]); #total number of cycles
                               63

These results agree with VV's explicit list of the cycles of this graph.

My little procedure above will gain a huge efficiency improvement by using remember tables for the traces, but I don't have time right now. If someone else chooses to undertake this, recall that a remember table cannot be indexed correctly by a mutable structure like a Matrix. So, use S and k as the table indices.

Here's a procedure a little more general (because it works for an arbitrary number of vectors) than what you asked for:

InSpan:= proc(V::{{set,list}(Vector[column]), Vector[column], Matrix}, w::Vector[column])
   try LinearAlgebra:-LinearSolve(`if`(V::Matrix, V, `<|>`(`if`(V::Vector, V, V[]))), w)
   catch "inconsistent": return false
   end try;
   true
end proc:

So, the first argument can be a column vector, a matrix, or a set or list of column vectors; and the second argument is a column vector. Example:

InSpan({ <1,2>, <3,4> }, <5,6>);

Regarding the terminology: The span of a set of vectors is the vector space of all linear combinations of those vectors. 

Several points to mention:

1. The command with PDE(tools) is nonsense to Maple. You mean with(PDEtools);. However, I encourage to not use this unless we determine that it's necessary. (In particular, PDEtools can be used to obscure the functional dependencies, which might be nice for the display of these equations, but makes it impossible for me to understand them and help you further.)
2. Using unprotect on a name doesn't make it safe to use that name! Change unprotect(gamma); to local gamma;.
3. Square brackets must not be used in place of parentheses.
4. I need to see the functional dependencies is order to help you further. For example, if v__1 is a function of x and t then you must write v__1(x,t) wherever they are used (or you can use v__1(t,x), as long as you use the same variable order consistently throughout the problem.

Your revised non-recursive procedure takes a long time due to expression swell, which is the tendency of expressions to get larger (literally larger, in terms of memory required, not numerically larger) when iterative operations are performed with exact arithmetic. If you change the first argument of your procedure from 10 to 2, you'll see what I mean. Newton's method in particular will usually cause an extreme amount of swell (in exact arithmetic). The solution is simple: Use floating-point (i.e., decimal) arithmetic instead. You can do this by wrapping every computation with evalf(...) (short for evaluate in floating-point arithmetic). So, here's a revised procedure:

Newton:= proc (n, y0, Q, b, m, k, S0) 
local A, P, Sf, y, F, dAdy, dPdy, Fdiff, i;
   y[0] := evalf(y0);
   for i from 0 to n do
      A[i]:= evalf((b+m*y[i])*y[i]);
      P[i]:= evalf(b+2*y[i]*sqrt(1+m^2));
      Sf[i]:= evalf(Q*abs(Q)*P[i]^(4/3)/(k^2*A[i]^(10/3)));
      F[i]:= evalf(S0/Sf[i]-1);
      dAdy[i]:= evalf(b+2*m*y[1]);
      dPdy[i]:= evalf(2*sqrt(1+m^2));
      Fdiff[i]:= evalf((2/3)*(F[i]-1)(5*dAdy[i]/A[i]-2*dPdy[i]/P[i]));
      y[i+1]:= evalf(y[i]-F[i]/Fdiff[i])
   end do;
   [seq(y[i], i = 0 .. n)]
end proc;

I put evalf in every assignment statement. Usually just a few (or even one) well-placed evalfs are sufficient, but the extra-time used by having unnecessary evalfs is too small to measure.

I notice that y[1] occurs inside your loop. I think that that should be y[i].

There is a residue command. If we take for granted that the enclosed poles for both integrals are {-I, I}, then the second integral is

(2*Pi*I)*add(map2(residue, f(2,z), z=~ {I,-I}));

and the first integral is

(2*Pi*I)*add(map2(residue, f(sqrt(2),z), z=~ {I,-I}));

which agrees with the answer you got from the int command.

It would be nice if there was an easy and general way to know that those are the relevant poles, but nothing comes immediately to mind. The singular command easily finds all the poles, but I don't know a trick to determine whether they are enclosed by the contour.  (But an idea is formulating in my mind.)

You need to prune off the spurious imaginary parts as soon as they appear, not after the computation is done. If not, the trajectory is like that of a spaceship that's off by a small fraction of a degree: If not corrected, it could miss its target by millions of kilometers.

Just curious: Why do you choose the 5th fixed point out of the 32 fixed points (2^5 = 32) of the 5th composition of the map?

Here's some modified code:

restart:
Digits:= trunc(evalhf(Digits)):
NestList:= proc(f, x, n::nonnegint)
local k, R:= hfarray(0..n, [evalhf(Re(x))]); 
   evalhf(
      proc(f, R, n) local k; for k to n do R[k]:= evalhf(f(R[k-1])) od end proc
      (f, R, n)
   ); 
   R
end proc:
logistic:= x-> 4*x*(1-x):
fp:= [solve((logistic@@5)(x) = x, x)];
cfp:= allvalues(fp[5]); #Why 5 _here_?
nstep:= 12:
p:= nstep+50; 
aaa:= NestList(logistic, cfp, p); 
dat:= <<[$-nstep-1..49]> | aaa^%T >: 
plot(
   dat, 
   style= pointline, symbol= solidcircle, symbolsize= 4, 
   thickness= 0
);

You can see that now the pattern holds for 9-10 periods, then starts to breakdown due to floating-point considerations.

By the way, I hate seeing my beautiful, painstakenly hand-formatted code converted into a "brick" of crap by being copy-and-pasted in 2D Input and then copy-and-pasted back to plaintext.

In Maple, it's a one-liner, requiring no local variables:

ProxySort:= (S::list,P::list)-> S[sort(P, 'output'= 'permutation')]:

The expression sqrt(x)/x is unstable in Maple. As soon as it's  created (whether by directly typing it or by rationalize), it gets automatically simplified to 1/sqrt(x). 

Like this:

MyMatrix:= apply~(Matrix((2,3), symbol= elem), t);

I cannot get LinearFit to fit a single constant parameter. It seems that that case is explicitly forbidden. But it can be easily done with plain minimize:

f:= x-> 3*x+a:
S:= minimize(add((f~(pts1)-pts2)^~2), location);
plot(
   [<pts1|pts2>, eval(f(x), indets([S], identical(a)=anything))], 
   x= (min..max)(pts), 
   style= [point,line], symbolsize= 20
);

Which gives a = -61/3, same as Kitonum's method.

Use

simplify(evalc([(()->"Re")=Re, (()->"Im")=Im](Ec))) assuming positive;

Hint: Substitute epsilon= 1/upsilon, then simplify, then consider "by hand" the limit as upsilon approaches +infinity. I think that this is easier to see than the limit as epsilon approaches 0 from the right. Likewise, limit(..., epsilon= 0, left) is equivalent to limit(..., upsilon= -infinity).

And, of course, to make your problem more understandable to Maple, you must replace [ ] with ( ) and 0.5 with 1/2 (as VV said).

Here is an extendable template that should let you answer (in a natural way if you know statistical language) any such question for any discrete distribution with integer support. All results come directly from applying simplify (with assuming) and sum to the pmfs of the distribution. Each distribution is represented by a Maple container structure called a Record, which contains fields for the pmf, support, and assumptions on the parameters.

All parameters can be specified as any algebraic expressions (including constants); they need not be simple names. The `if` expressions in the Records adjust the assumptions appropriately.

The table (actually a DataFrame) at the end actually looks like a table when viewed in a worksheet. The MaplePrimes renderer squashes it here.
 

Download DiscreteDist.mw