In my first response I thought the question was the same as the OP asked me privately.
This response is now modified some.
Notice that I made the variables F0,F1,F2 local. The secant metod (and other methods) evaluates p at various points pp2. It is not necessarily the last computed p(pp2) that gives the desired pp2. However, if F0,F1,F2 are global then they will correspond to the last pp2 used.
Also be careful about assigning to p2 if you are running the procedure p again: it contains a subs(p2=pp2, ...).
Actually this might in this case be the real problem.
The answer is to use locals and not to assign to p2.
p:=proc(pp2) local res,F0,F1,F2; etc. end proc:

Putting a print statement into the solution procedure when using the procedural version of dsolve/numeric/IVP may be illuminating:

restart;
#An example
ode:=diff(x(t),t)=f(t,x(t));
f:=(t,x)->t+x^2;
#Solving the following problem:
ode:=diff(x(t),t)=f(t,x(t)); #with x(0)=1
dsolve({ode,x(0)=1});
#Now numerically
p:=proc(N,t,x,xp) print(t,x); xp[1]:=f(t,x[1]) end proc;
sol:=dsolve(numeric,number=1,procedure=p,procvars=[x(t)],initial=Array([1]),start=0,method=rkf45);
sol(.5);
#Compare with
sol:=dsolve(numeric,number=1,procedure=p,procvars=[x(t)],initial=Array([1]),start=0,method=classical[rk4],,stepsize=.05);

Doesn't unapply do what you want:

TestDerivative:=unapply(eval(diff(TestDerivative(y),y),y=x),x);

First of all there is a syntax problem: dsolve needs a set of equations as in dsolve({eq1,bcs1}).

You won't have any luck finding an analytic expression for your solution, so try numerically.
Here 'infinity' has to be replaced by a large enough number N:

N:=10:
bcs := (D(f))(0) = 1, f(0) = 0, (D(f))(N) = 0;
sol := dsolve({bcs, eq1},numeric);
plots:-odeplot(sol,[[eta,f(eta)],[eta,D(f)(eta)]],0..N);

With a little less work you can produce the same picture like this:
restart;
h := x->piecewise(x < 0, undefined, x <= 2, 5-x^2, 3-x):
plot([h(x),h(-x)], x = -10 .. 10, color = [red,green], thickness = 3);

Why did you set  z = 2?
By T(N^2(z)) do you mean T(N(z)^2) or do you mean T(N(z))^2?
In the latter case it is easy to see the general pattern since a^k - b^k has a-b as a factor for any positive integer k. In the former case it is easily seen that T(N(z)^2)-(T(z))^4 is not zero (even for z=2).
See the following:

restart;
T:=z->(z-I)/(z+I);
N:=z->(z^2-1)/(2*z);
simplify(T(N(z))-T(z)^2);
a:=T(N(z)):
b:=T(z)^2:
a1:=normal(a); #Rewriting a
b1:=normal(b,expanded); #Rewriting b
simplify(a1-b1);
simplify(a1^k-b1^k) assuming k::posint;
#Now the first interpretation:
T(N(z)^2)-(T(z))^4;
simplify(%); #Almost never zero
eval(%,z=2); #Not zero

With sys being the list of your 16 equations and noticing that the variable Ll (not L1) is missing from the variables set, you can use fsolve on sys[8,16] as follows:

g:={Lfd =0.233, R1d=0.0848, L1d=0.2, L1q=-0.09, Ladssec=0.099, Laqssec=-0.1, Rfd=0.00173, R1q=0.0129,Ll};
nops(g);
fsolve(sys[8..16],g);

returning
{L1d = .1989149191, L1q = -0.8847981237e-1, Ladssec = 0.9877218169e-1, Laqssec = -0.9817047856e-1, Lfd = .2305604744, Ll = 0.9815030978e-1, R1d = 0.8441614099e-1, R1q = 0.1292859817e-1, Rfd = 0.1715066697e-2}

You cannot use square brackets i.e. [ and ] for grouping. You must use parentheses i.e.  ( and ).
After making these changes the numeric solution as you have it works.

Your second question:

Try
res:=pdsolve({dotW11,icond});
plot(eval(rhs(res),t=1),y=22..26);
plots:-animate(plot,[rhs(res),y=22..26],t=0..1);

Several multiplication signs are missing, but that is not really the problem.
You can only collect with respect to names or unevaluated functions.
Also it is not clear to me what you want as outcome.
To find out where multiplication signs are missing you can do
indets(g,function);
#If that returns anything but the empty set you have forgotten a `*`.
You may try this
g := -2-k[1]*(lambda*alpha[2]*k[2]+alpha[1]*(-2-2*k[2]+k[2]*lambda^2))+k[1]*(lambda*alpha[2]*k[2]+alpha[1]*(-2-2*k[2]+k[2]*lambda^2))*lambda*(lambda*alpha[2]*k[2]+alpha[1]*(-2-2*k[2]+k[2]*lambda^2))^2;
indets(g,function); #Returns {}
u:=lambda*alpha[2]*k[2]+alpha[1]*(-2-2*k[2]+k[2]*lambda^2);
subs(u=w,g);
#And now the question is what do you really want to do?

If you write your system on the form  T'(t) = A.T(t) + b, where T = <T[1], ..., T[n] > and b is  a constant vector, then you can find the eigenvalues  of A easily by LinearAlgebra:-Eigenvalues.
I tried with n = 20 and found 19 negative and 1 positive eigenvalue. The latter is likely to create problems for larger times even if you start with initial conditions in the space spanned by the eigenvectors corresponding to the negative eigenvalues because of roundoff errors.
I found the same problem for the two other values I tried n=2 and n = 5.

To avoid retyping you can do:
Lsys:=convert(sys,list);
A1,b:=LinearAlgebra:-GenerateMatrix(subs(seq(diff(T[i](t),t)=0,i=1..n),Lsys),[seq(T[i](t),i=1..n)]);
A:=-A1;
A.<seq(T[i](t),i=1..n)>+b; #Vector of right hand sides in T'(t) = A.T(t) + b
LinearAlgebra:-Eigenvalues(A);

Added:
I tried not assigning to your constants except n.
Then I did:
LinearAlgebra:-IsDefinite(A,query='negative_definite');
simplify(%) assuming positive;
and examined these inequalities with the specific values you have chosen for the constants.
The problem seemed to be the last inequality. The choice of lambda is essential.

You need a  list of expressions when plotting several expressions:

plot( [exp(x),exp(x)+3,exp(x)*2,exp(x+3)] ); #Maple finds the domain itself
plot([exp(x),exp(x)+3,exp(x)*2,exp(x+3)],x=-6..3.5); #You give the domain
#No list is  required for one expression
plot(ln(x),x=0..15);

Look back to the answer I gave to your question at
http://www.mapleprimes.com/questions/152986-How-To-Put-The-Condition-In-While-z

I reprint the example in one of my comments:
restart;
Orbit:=proc(Map,ic,Nmax::posint:=100)
  local orbit,z,t;
  orbit:=Vector();
  orbit(1):=ic;
  z:=ic;
  for t from 2 to Nmax while z<>ic or t=2 do
    z:=Map(z);
    orbit(t):=z
  end do;
  orbit
end proc;
f:=L->modp~(L+[8,9],5);
orb:=Orbit(f, [1, 1], 15);
plot(orb,style=point,symbol=solidcircle,symbolsize=20);

M:=Matrix([[1,2,3,4,5],[6,7,8,9,10]]);
L:=convert(M,listlist);
(cat@op)~(L); #Now symbols
convert~(%,string);

Try this instead.

restart;
with(plots):
MW := 29; #assignment uses :=
cp := 1200;
hf := 4*10^7;
AF := 16;
eq1 := (D(F))(t) = 6.19*10^9*exp(-15098/T(t))*F(t)^.1*(.2*O(t))^1.65; #T(t)!!!! not T
eq2 := (D(O))(t) = 16*(D(F))(t);
eq3 := (D(Pr))(t) = -17*(D(F))(t);
eq4 := (D(T))(t) = (D(F))(t)*hf*(8314*T(t))/((cp-8314)*P(t));
eq5 := (D(P))(t) = P(t)*(D(T))(t)/T(t)-P(0)*(D(T))(t)/T(0);
ics:={O(0) = 152.32*10^(-3), P(0) = 101325, T(0) = 753, F(0) = 9.52*10^(-3), Pr(0) = 0};
#Equations eq2 and eq3 ought to be left out since O(t) and Pr(t) are trivially found from F(t):
#eq2:
eval(O(0)=16*F(0)+C,ics);
solve(%,C);
#eq3
eval(Pr(0)=-17*F(0)+C,ics);
c:=solve(%,C);
#Thus
OF:=O(t)=16*F(t);
PrF:=Pr(t)=-17*F(t)+c;
#The reduced system:
ics2:={P(0) = 101325, T(0) = 753, F(0) = 9.52*10^(-3)};
sys2:=eval({eq1, eq4, eq5},{OF,PrF} union ics2);
sol := dsolve(sys2 union ics2, numeric);
odeplot(sol,[[t,P(t)]],0..10^6);
odeplot(sol,[[t,F(t)]],0..10^6);
odeplot(sol,[[t,T(t)]],0..10^6);

A:=Matrix([[8,4,2,1],[4,8,2,1],[2,2,8,1],[1,1,1,8]]);
with(LinearAlgebra):
IsDefinite(A);
L,U:=LUDecomposition(A,output=['L','U']);
D1:=DiagonalMatrix([seq(U[i,i],i=1..4)]);
L.D1.Transpose(L);