It might help if you tried to explain what you are trying to do.
Your a and b are both found to be the 3 real numbers 0., .9865070072, 1.013677544. You are asking complexplot to plot {0., .9865070072, 1.013677544,0., .9865070072, 1.013677544}.
In fact it is doing that: Try adding thickness=3 to the command. complexplot joins the points with straight lines. Compare with complexplot({0,1+I,1-I});
Notice also that because you use curly brackets (for a set) your set
{0., .9865070072, 1.013677544,0., .9865070072, 1.013677544}
is first reduced to
{0., .9865070072, 1.013677544}.
Your use of complexplot3d before the command with(plots) obviously doesn't produce anything, but after with(plots) it does.

Begin by changing all occurrences of pure names A, B, phi into the functional form A(r), B(r), phi(r).
To make sure you have done it all over test with
indets(odesys,name);

until none of A,B,phi show up.

If you call the result coming out of odetest for res then try
res1:=eval(res,{a=1,alpha=1,beta=2,gamma=3,lambda=4,w=1});

and then

simplify(res1);

You will see that you don't get an expression which is identically zero.

If you let y1 denote the result of the last subs you make, then try:

plots:-complexplot3d(x*y1,x=-2-I..2+I);
fsolve(x*y1,x);
fsolve(x*y1,x=1);

Take a look at timelimit.

?timelimit

When I run your code in Maple 2015 I don't get any complaints about singularity.
Be aware that fieldplot3d and display belong to the plots package.

However, if you extend the t-interval to the left e.g. to -6 then you get warnings about singularities.

You could try the first initial condition in IC and use dsolve/numeric like the following.
The construction {SYS[],IC[1][]} just gives you a set containing the odes and the initial condition.

res:=dsolve({SYS[],IC[1][]},numeric);
plots:-odeplot(res,[t,y(t)],-6..5);
res(-5.1383855);

#There is most likely indeed a singularity at about t=-5.1383855. What this means here is that the orbit "reaches" infinity in a finite time. There is nothing weird about that.
The simple example dsolve({diff(y(t),t) = -y(t)^2, y(0)=1}) has a singularity at t=-1.

I tried putting m0= 10*n0 and then animate with n0 = 1.. 25. For the rest of the parameters I used your values.
What appears from the animations is that for n0 < 4.633... the point (0,0) is the only equilibrium point and it is an attractor (a node).
For 4.633 < n0 < 15.699 two more equilibria appear and (0,0) is a saddle point.
For n0 > 15.699 an additional two equilibria appear and (0,0) is now a repelling node.
I believe that no qualitative changes take place after n0 > 15.699 and that any problems after that are numerical. Your n0-value is 1000, much larger than 25 of course.

I only drew one curve and for the sake of visibility I took x(0)=0, y(0)=0.5 instead of your initial condition.
implicitplot is used to draw the nullclines.

restart;
Digits:=15:
   u := am+bm*m0*x(t)+cm*n0*y(t);
   v := an+bn*m0*x(t)+cn*n0*y(t);
   eq1 := diff(x(t), t) = 2*A*(sinh(u)-x(t)*cosh(u));
   eq2 := diff(y(t), t) = 2*A*(sinh(v)-y(t)*cosh(v));
#param changed:  
param := A = 1/2, am = 0, bm = 1.2*exp(-4), cm = .5*exp(-5), m0 = 10*n0, an = 0, bn = -exp(-4), cn = 1.2*exp(-3);
#Animating nullclines:
RHS:=evalf(subs(param,x(t)=x,y(t)=y,rhs~([eq1,eq2])));
plots:-animate(plots:-implicitplot,
    [RHS,x=-2..2,y=-2..2,gridrefine=3,thickness=2,color=[red,blue]],
    n0=1..25,frames=50
    );
p1:=%:
#Animating phaseportrait with one curve:
plots:-animate(DEtools[DEplot],
    [subs(param,[eq1,eq2]),[x(t),y(t)],0..10,[[x(0) = 0, y(0) = 0.5]],
      x=-2..2,y=-2..2,color=green,arrows=small,linecolor=black,thickness=3,stepsize=.01],
     n0=1..25,frames=50
     );
p2:=%:
#Displaying the two animations simultaneously:
plots:-display(p1,p2,size=[1000,600]);
#Making the jacobian a function of n0 , x, y:
J:=unapply(VectorCalculus:-Jacobian(RHS,[x,y]),n0,x,y):
#Check that (0,0) is an equilibrium point for all values of n0:
eval(RHS,{x=0,y=0});
#Determining the type and stability of the point (0,0):
LinearAlgebra:-Eigenvalues(J(n0,0,0));
LinearAlgebra:-Determinant(J(n0,0,0));
solve(%=0,n0);
#Trying n0=5 where there are 3 equilibria. Finding one of these:
fsolve(eval(RHS,n0=5),{x,y});
#Determining the type and stability of that point:
J(5,op(subs(%,[x,y])));
LinearAlgebra:-Eigenvalues(%);
#Trying n0=17 where there are 5 equilibria. Finding one of these:
fsolve(eval(RHS,n0=17),{x=1,y=-1});
#Determining the type and stability of that point:
J(17,op(subs(%,[x,y])));
LinearAlgebra:-Eigenvalues(%);

#There is obviously more to be done, but this must suffice for now.

The response from Maple to your code (ignoring the save statements, the attempt at at exact solution of the ode, and the first assignment to init) is:
Warning, cannot evaluate the solution further right of .17412248e-1, probably a singularity

That seems indeed to be the case: (see note below)

Try
odeplot(sol, [[t,x(t)],[t, y(t)]], 0 .. 0.017412248);
You will notice that the solution curve for x(t) shoots "straight up" and that for y(t) "straight down" at (x,y) values given by  the result of

sol(0.017412248);

This is just how it is.

Edit: I was too hasty in my reply. As can be seen from Rouben Rostamian's answer the problem is a numerical one.

You cannot use list brackets, i.e.  [ ]   instead of parentheses. So remove those and replace them with ( ).

I don't know what you are referring to with dsnumsort, C1, C2, but you can plot the orbit and the time dependence like this:

restart;
dx := diff(x(t), t, t) = -G*Mz*x(t)/(x(t)^2+y(t)^2)^(3/2);
dy := diff(y(t), t, t) = -G*Mz*y(t)/(x(t)^2+y(t)^2)^(3/2);
G := 6.67*10^(-11); Mz := 6*10^24:
IniC := x(0) = 7*10^6, (D(x))(0) = 0, y(0) = 0, (D(y))(0) = 9*10^3;
Digits := 15:
Ns := dsolve({dx, dy, IniC},numeric);
 
plots:-odeplot(Ns,[[t,x(t)],[t,y(t)]],0..16000); #Time dependence
plots:-odeplot(Ns,[x(t),y(t)],0..16000,scaling=constrained); #Orbit
plots:-odeplot(Ns,[x(t),y(t)],0..13060,scaling=constrained,frames=100); #Orbit animated
#You can even make dsolve stop when one period has been finished:
NsE := dsolve({dx, dy, IniC},numeric,events=[[[y(t)=0,x(t)>0],halt]]);
plots:-odeplot(NsE,[x(t),y(t)],0..16000,scaling=constrained);

Your problem is that your boundary conditions involve a derivative that is not normal to the boundary. Although it is somewhat hidden by your use of v instead of D[2](R) it is a condition such as
a*v(-30,t) = b* D[1](R)(-30,t)
which is the problem.
In reality it is
a*D[2](R)(-30,t) = b* D[1](R)(-30,t)
This condition involves the derivative a*D[2](R)  at the boundary l = -30.
This derivative is not normal to the boundary. (See error message at the bottom of this).

Try this.
Define as suggested by Carl Love (and leaving out the last condition on u)

bdry00 := ((30^2+1^2)/30^2)^(1/4)*v(-30, t) = -((30^2+1^2)/30^2)^(1/2)*(D[1](R))(-30, t), ((30^2+1^2)/30^2)^(1/4)*v(30, t) = -((30^2+1^2)/30^2)^(1/2)*(D[1](R))(30, t), u(-30, t) = 0;

Assume that PDE01 and PDE03 have been defined as you did. Then do:

PDE01a:=subs(v(l,t)=diff(R(l,t),t),PDE01); #Using your PDE02
bdry:=simplify(subs(v=D[2](R),{bdry00})); #Rewriting bdry00
#Also rewriting your init00:
inits:={D[2](R)(l, 0) = 0, R(l, 0) = 0, u(l, 0) = sqrt((l^2+1^2)^(1/2))*10^(-5)*exp(-(l-10)^2/.5^2)};
pdsolve({PDE01a,PDE03}, bdry union inits, numeric);
## You get the error message:

Error, (in pdsolve/numeric/process_IBCs) initial/boundary conditions can only contain derivatives which are normal to the boundary, got (D[2](R))(-30, t)

I don't know what exactly you are comparing, but I notice that you are missing the term 1 of s in the unnamed sum.
If you remember that the difference between your two ways of calculating s at s=70 is -3.*10^(-1394)

restart;
Digits:=1400:
p:=convert(series(cos(x),x,201),polynom):
s:=add(op(i,p),i=[1,2,3,-1,-2,-3]);
##Here I added the missing 1:
r1:=1+add(evalf(coeff(p,x,i))*70^i,i=[200,198,196])+add(evalf(coeff(p,x,i))*70^i,i=[2])+add(evalf(coeff(p,x,i))*70^i,i=[4]);
simplify(eval(s,x=70));
r:=evalf(%);
r-r1;

Notice that what you write down as (the) solution is actually a sequence of three "solutions".
The reason for my use of quotes is that it is questionable whether the last two are solutions. See bottom.

The ode, s, is autonomous so we may as well assume that the initial condition is given at t=0.

The "solutions" all satisfy h(0)=14/5 (=2.8).  That is a rather special initial condition, which makes the right hand side of the ode undefined at t=0. So expect problems.

The same 3 solutions you get by finding the general solution by dsolve(s); and then setting _C1 = 0 after which you solve for h(t).
Call the sequence of solutions sol.

The first solution sol[1] is in fact a solution on the interval 0 .. infinity.
odetest(sol[1],s); #returns zero, meaning that it is indeed a solution
plot(rhs(sol[1]),t=0..25); #Notice the slope at t=0.

#If you insist on plotting all 3 in the complex plane:

plots:-complexplot([seq(rhs(sol[i]),i=1..3)],t=0..25,thickness=3);

Notice that sol[2] and sol[3] are imaginary for t=0..25.

More illuminating is the animation:
plots:-animate(plots:-complexplot,[[seq(rhs(sol[i]),i=1..3)],t=0..T,thickness=3],T=0..25);

## Finally we can see that the two last "solutions" are not solutions:
res23:=odetest~([sol[2..3]],s); #Doesn't immediately (at lest) return [0,0].
#Let us animate:
plots:-animate(plots:-complexplot,[res23,t=0..T,thickness=3],T=0..25);
#Not at all zero on any part of the interval t=0..25

Maybe you could use events as in this very simple example, where events=[[y(x)=1e-10,y(x)=0]] means that if y(x) during integration attains the value 1e-10 then the action taken is that y(x) is set to zero, which in this case means that it remains zero.

ode:=diff(y(x),x)=-10*y(x);
sol:=dsolve({ode,y(0)=1},numeric,events=[[y(x)=1e-10,y(x)=0]]);
sol(4);
plots:-odeplot(sol,[x,y(x)],0..4);

## without events:
sol2:=dsolve({ode,y(0)=1},numeric);
sol2(4);

########  However, consider this example:
ode2:=diff(y(x),x)=-10*y(x)+piecewise(x<4,0,1);
sol3:=dsolve({ode2,y(0)=1},numeric,events=[[y(x)=1e-10,y(x)=0]]);
sol3(5);
plots:-odeplot(sol3,[x,y(x)],0..5);

## So maybe you will want the integration to stop instead:
sol4:=dsolve({ode2,y(0)=1},numeric,events=[[y(x)=1e-10,halt]]);
sol4(4);

You have a syntax problem:  B.C should be BC.

restart;
eq := diff(T(r), r, r)+(diff(T(r), r))/r+(3*n+1)*Nu*(1-r^((n+1)/n))*T(r)/(n+1) = 0;
dsolve(eval(eq, n = 1));
BC := D(T)(0) = 0;
dsolve(eval({eq,BC},n=1));
## Kummer behavior near zero:
plot(eval([KummerM(1/2-(1/4)*sqrt(2)*sqrt(Nu), 1, sqrt(2)*sqrt(Nu)*r^2),KummerU(1/2-(1/4)*sqrt(2)*sqrt(Nu), 1, sqrt(2)*sqrt(Nu)*r^2)],Nu=1),r=0..1);

#######
The ode has a singularity at 0.

The general solution is found by dsolve as a linear combination of KummerM and KummerU multiplied by exponentials.
The second term tends to infinity as r - > 0 (right). Its derivative tends to -infinity.
The first term tends to 1 as r->0. Its first derivative tends to 0.
The conclusion is that _C2 has to be zero if you want D(T)(0)=0 and _C1 can be chosen arbitrarily.

You are using numerical solution with the syntax you are using. Thus s needs to have a concrete value.

The output is a module exporting plot, plot3d, animate, value.

See the help page
?pdsolve,numeric

Example:
pds := pdsolve(eval(PDE,s=1), IBC, numeric);
pds:-animate(z=3,frames=50);