Why do you say the answer is wrong?  What do you expect?  0?  If so, don't ask for a series but a limit:

limit(f, x=infinity);
                                0

For the third task, it is more efficient to use PolynomialTools:-CoefficientList (or CoefficientVector):

PolynomialTools:-CoefficientList(convert(A,polynom),x);      
                                                                 -1
                                   [0, 1, 0, -1/6, 0, 1/120, 0, ----, 0, 1/362880]
                                                                5040

Do you mean why the word "unapply" and not some other name?  The reason, I believe, that "unapply" was chosen is that the operation is, in some sense, the inverse of function application (which can be done with the "apply" command).  See the help page for apply, which states

- Using unapply, the apply procedure satisfies the functional equation 
  apply(unapply(p, v), op(v)) = p. 

Can you give an example?  Have you tried executing the examples in the BodePlot help page (use the button that converts a help page to a worksheet)?

I wouldn't suggest the following---it presumes some stuff that I'll guess you haven't learned---but others may find it interesting:

T := 4*x^2 - 4*x*y + y^2:
C := x^2 + y^2 - 5^2:

with(difforms):
defform(x=scalar, y=scalar):
dTC := factor(scalarpart(d(T) &^ d(C)));
                         dTC := 4 (x + 2 y) (2 x - y)

sol := solve({dTC,C}, [x,y]):
map(subs, sol, T):
convert(%,radical);
                                   [125, 0]

Nicer, but also not suitable, I suspect, is

with(Student[MultivariateCalculus]):
stationary := LagrangeMultipliers(T,[C],[x,y],output=detailed):
convert(map(subs, [stationary], T), radical);
                                   [0, 125]]

Its not clear what you want.  Just a plot?  Or are you trying to fit curves to the data?

While it may not solve your problem, you can compute the asymptote of the transformed integral:

evalindets(eq3, Int(anything$2), limit, k=infinity);
                                      1
                                      ---
                                      4 k

If all you are interested in is the assignments, then it should be clear

(1) 2
(2) n+1 [i is assigned n time] 
(3) n [once each iteration, so n total]
(4) n [once each iteration, so n total]
(5) 0
(6) 0
So the total is 2*n+3.  The reason that (2) is n+1 and not n is that on the last iteration i is assigned to n+1, which terminates the loop.

The complexity of this algorithm, however, is not necessarily O(n), at least not in Maple.  It depends on what y is assigned.  Assume it is a not assigned.  To see what happens do

printlevel := 20:
n := 3:
a:=1: b:=1:
for i to n do
    a := (a*y)/i;
    b :=b+a;
end do;
                                    a := y

                                  b := 1 + y

                                          2
                                         y
                                   a := ----
                                         2

                                                2
                              b := 1 + y + 1/2 y

                                          3
                                         y
                                   a := ----
                                         6

                                           2        3
                         b := 1 + y + 1/2 y  + 1/6 y

Note that b increase in size with each loop.   To better see this you can use the length function as a crude measure of the complexity of intermediate expressions:

printlevel := 1:
data := table():
n := 1000:
a:=1: b:=1:
for i to n do
    a := (a*y)/i;
    b :=b+a;
    data[i] := [length(a), length(b)];
end do:

with(plots):
display(loglogplot([seq([i,data[i][1]], i=1..n)]),
        loglogplot([seq([i,data[i][2]], i=1..n)]));

If you plot this you will see that the length of a increases linearly with n, while the length of b increases as the square of n.  Because all intermediate values are saved in the Maple simplification table, the complexity of the entire algorithm could be O(n^3).  Now rerun this test, but assign y a numeric value, say 1.  You will then see that both a and b increase linearly.  Finally, rerun once more, but with y assigned a float, say 1.0.  Now the complexity (as measured by length) of a and b does not increase with n.

Robert's use of simplifying with a  side relation, see ?simplify/siderels, is elegant and practical.  If you didn't know about that, you could, at least for these simple problems, solve for x, pick one solution, and then substitute into the result:

sol := solve((x+1)^2 = a, {x}):
normal(subs(sol[1], (x^2+2*x+1)^3+2)); 
                                         3
                                    2 + a

Edited

That same method, but using expand, will work for your followup example (exp is a function in Maple, not a symbol):

y := exp(3*x+3)+2; 
sol := solve(exp(x+1)=a,{x});
expand(subs(sol,y));      
                                         3
                                    2 + a

The command display is part of the plots package.  Access it either with the long name, plots:-display, or execute with(plots) before calling display.

The assignment

y1 := x -> y1(x);

and the subsequent call to y1 is the culprit.  The easy way to figure this out is to do

debug(PS):
PS();
.... [output omitted]
                               y1 := x -> y1(x)

<-- ERROR in PS (now at top level) = "too many levels of recursion"}
Error, (in y1) too many levels of recursion

What you should do is

y1 := 'y1';

Similarly for y2, later in the procedure. Finally, using with(plots) inside a procedure does not work.  What you want to do is remove that statement and change display to

plots:-display(...);

With those changes I get a plot.

This might get you started,

q := sqrt(2*m*(V-x))/h;
z := ( exp(-I*q*b)*(cosh(q*b)))^2 - I*sinh(q*b) * exp(-4*I*q*b);
simplify(evalc(abs(z)^2)) assuming V < x, m>0, h>0;

u := c*(exp(x-4*t)+exp(4*t-x))^(-2):
eq := diff(u,t) + u*diff(u,x) + diff(u,x$3) = 0:
sol := solve(eq,[c]):
simplify(sol);

Use dsolve/numeric with the output=array option, where array is an Array specifying the values of the dependent variable.  See ?dsolve/numeric. For example:

deq := {diff(y(x),x,x)+y(x)=0, y(0)=1, D(y)(0)=1}: 
dsolve(deq,numeric,output=Array([seq(0..1, 0.1)]));  

You've got the right idea, however, there are a few problems.  First, as assigned, d is merely an expression, not a procedure.  I'd write the procedure as

d := m -> numtheory:-sigma(m)-m; 

Next, the assignment to Amicable has a syntax error; the call to print should look like print(m,n);  Also, you don't want to make m a formal parameter.  What you probably should do is let the max number of loops be the formal parameter:

Amicable := proc(nmax)
local m,n;
   for m from 2 to nmax do
       n := d(m);
       if d(n) = m then
           print(m,n);
       end if;
   end do;
end proc:
Amicable(10000);

You might also want to add a test for m <> n, otherwise m=n, which isn't really a pair, but is a perfect number.
 

Note that this procedure merely prints the number.  More useful is to return the pairs in a list.  There are many ways to do that, but the easiest might be

[seq(`if`(d(d(n))=n and d(n)<>n, [n,d(n)], NULL), n = 1..10\000)];

A useful trick is to return the pairs as sets rather than lists, and as one big set (of sets), rather than a list, so that duplicates are automatically eliminated.