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

Joe has pointed out, essentially, that's it's often erroneous and usually suspicious to refer to a for loop's index variable outside the loop. The practice does have a few valid uses, such as determining whether the loop ended due to exceeding the maximum index or due to some other reason.

You have another major error: "Product of the eigenvalues" does not mean the dot product of the vector of eigenvalues with itself; rather, it means the (scalar) product of the three eigenvalues. You can do this with command mul, like this:

mul(N[i])

This is something like what your professor has in mind:

with(LinearAlgebra):
N:= 8:  n:= 3:  r:= true:
to N do
   M:= RandomMatrix(n,n);
   if Determinant(M) <> simplify(expand(mul(Eigenvalues(M)))) then 
      r:= false;
      break
   end if
end do: 
r;

I prefer this:

macro(LA=LinearAlgebra):
N:= 8:  n:= 3:  r:= true:
to N do
   if (radnormal@(LA:-Determinant - mul@LA:-Eigenvalues))(LA:-RandomMatrix((n,n))) <> 0 then
      r:= false;
      break
   fi
od:
r; 

I prefer it primarily because I trust radnormal when applied to a zero-equivalent radical expression more than I trust simplify@expand when applied to an integer-equivalent radical expression.

If F is list (or set or sequence) of functions (or procedures) and S is a sequence of arguments that would work for any of them, then simply F(S) is what you want. Note that this uses neither map nor an elementwise operator with ~. Example:

[iquo, irem](7,5)

or

(iquo, irem)(7,5)

I often use this as (min, max)(S). So, the number of arguments doesn't matter.  In particular, it could be one argument.

If X and Y are lists of the same length, where X are arguments for the first position and Y for the second, then do

F~(X, Y)

Unlike zip, this idea can be extended to any number of arguments. Any of the argument lists can be replaced by a single element, as long as that element is not itself a list, set, table, or rtable. The take-away lesson from all this is that there's no syntactic difference between a single function (or procedure) and a list, set, or sequence of them (an explicit sequence needs to be in parentheses).

Your code

map(eval~, [f(x), g(x)], x=~ [p,q,t])

can be (and should be) replaced by 

map(map, [f, g], [p, q, t])

because it's not necessarily possible to evaluate every function at symbolic x, and even if it's possible, it's not necessarily efficient.

The following procedure detects and reports cycles of exact algebraic numbers:

CycleDetectAlgnum:= proc(f, x0::radalgnum, Max::posint:= 999)
local x:= evala(Normal(x0)), R_inds:= table([x=0]), R_vals:= table([0=x]), k, j; 
   for k to Max do
      x:= evala(Normal(f(x)));
      if assigned(R_inds[x]) then return seq(R_vals[j], j= R_inds[x]..k-1) fi; 
      R_inds[x]:= k;  R_vals[k]:= x
   od;
   return
end proc
:   

It's usage is exactly like the other's:

CycleDetectAlgnum(y-> 4*y*(1-y), (5-sqrt(5))/8);

If you're only interested in detecting and reporting cycles rather than analyzing algebraically why they happen, then the following code---which does it all in hardware complex floats---will be much faster than either iterating the map or using a closed form from rsolve.

CycleDetectInner:= proc(
   C::Array(datatype= complex[8], order= C_order),
   k::integer[4],
   x::complex[8],
   eps::float[8]
)::integer[4];
option autocompile;
local j::integer[4];
   for j to k do
      if abs((C[j]-x)/`if`(abs(x) <= eps, 1, x)) <= eps then return j fi
   od;
   0
end proc
:
CycleDetect:= proc(
   f::procedure, 
   x0::complexcons, 
   Max::posint:= 999, 
   eps::{float,positive}:= evalhf(DBL_EPSILON)
)
option hfloat;
local x:= evalhf(x0), C:= Array(1..1, [x], datatype= complex[8], order= C_order), k, p; 
   for k to Max do
      x:= f(x);
      p:= CycleDetectInner(C, k, x, eps);
      if p > 0 then return C[p..] fi;
      C(k+1):= x
   od;
   return
end proc
:   

Example usage:

CycleDetect(y-> 4*y*(1-y), (5-sqrt(5))/8);

You may want to apply Re~ to the results to remove the zero imaginary parts. Any returned result is a 1-D Array containing a detected cycle. A return of NULL means that Max iterations were performed without a repeated value. 

The reason that add doesn't work as expected is that Tolerances doesn't overload it, whereas it does overload `+`. So, if you replace add with (`+`@op), you'll get your expected results.

Another option is post-processing of the partially evaluated results returned by add:

subsindets(Z, :-`+`, `+`);

The second `+` in that command is properly Tolerances:-`+`, but under the auspices of with(Tolerances), simply `+`  is sufficient.