This is weird indeed, but has clearly to do with the appearance of a''(t) as a 'known' function.

There is a similar procedure for bvp-problems `dsolve/numeric/bvp/convertsys`.
When you try

`dsolve/numeric/bvp/convertsys`(test, [], x(t), t, y, yp );

you receive the error message

Error, (in dsolve/numeric/bvp/convertsys) specified dependent variables x(t) do not agree with input system {a, x}

A workaround is to spell diff incorrectly as in ddiff(a(t),t,t) or of course just use a notation like app(t).

Once you receive the result from convertsys, you can easily replace those by what they ought to be.
################
The question remains why replacing 'a' by 'p' helps.

Looking at the procedure `DEtools/convertsys` in particular at the loop starting with line 51:
for idiff in diffset do ....
it appears that the order in which the elements of diffset are taken is different in the two cases (a(t) versus p(t)) in Maple 2015 and also in Maple 18, but not in Maple 12 or 15.
I found this confirmed also by
diffset:={diff(a(t), t), diff(x(t), t), diff(a(t), t, t), diff(x(t), t, t)};
convert(diffset,list);
diffset:={diff(p(t), t), diff(x(t), t), diff(p(t), t, t), diff(x(t), t, t)};
convert(diffset,list);
In Maple 2015 (and 18) the a's come first in the first case. In the second case the second derivatives come first.
In Maple 12 and 15 the first derivatives come first and convertsys doesn't work for either use of letter.

You have 'res' in the definition of T0 and T1 inside the procedure Q2. That line should be:

T0,T1:=op(subs(subs(solT),[theta(x),diff(theta(x),x)]));

You may consider simplifying EQT2 by doing:

EQT2:=simplify(EQT2);

That makes EQT2 shorter.

You could do as follows:

Th0:=[seq(cat(Y1,i),i=1..5)](0); #The theta(0) values
PrTh0:=zip(`[]`,L,Th0); #The points in the Pr-theta(0) plane
plot(PrTh0,labels=[Pr,theta(0)]); #Points connected
plot(PrTh0,labels=[Pr,theta(0)],style=point,symbol=solidcircle,symbolsize=20); #Points not connected

#################################
Using an arbitrary number of Pr values:

restart;
  Digits:=15:
  Eq1 := diff(f(eta), eta, eta, eta)+f(eta)*(diff(f(eta), eta, eta))-(diff(f(eta), eta))^2 = 0;
  N := 1;
  blt := 10;
  Eq2 := (diff(theta(eta), eta, eta))/Pr+f(eta)*(diff(theta(eta), eta)) = 0;
  bcs1 := f(0) = 0, (D(f))(0) = 1, (D(f))(blt) = 0;
  bcs2 := (D(theta))(0) = -N*(1+theta(0)), theta(blt) = 0;
# Procedure to compute theta(0) for given value of Pr:

p:=proc(pr) local res;
  res:=dsolve(eval({Eq1, Eq2, bcs1, bcs2}, Pr = pr), numeric, output = listprocedure,maxmesh=1024,initmesh=1024);
  subs(res,theta(eta))(0)
end proc;
 
plot(p,2.5..10);

p:=x^(n-1)+y^g+x^m+5*x^(n-2);
max(op(map2(op,2,indets(p,`^`(identical(x),anything)))))  assuming n>=m+1;

It is sometimes better (in my view) to use linear algebra for linear systems. It provides more insight.

I use the version of your system as corrected by Markiyan.

restart;
sys_ode := diff(ph(t),t) = (1-yc)*pc(t)+yh*prj(t)+urd*prd(t)+ugd*pgd(t)-yc*ph(t), diff(pc(t),t) = yc*ph(t)-(2-yc)*pc(t), diff(pa(t),t) = ya*pc(t)-pf(t), diff(prj(t),t) = yrj*pa(t)+prj(t), diff(prd(t),t) = -urd*prd(t)+yrd*pa(t), diff(pgd(t),t) = -ugd*pgd(t)+ygd*pa(t), diff(pf(t),t) = (1-ygd-yrj-yrd)*pa(t)+(1-yh)*prj(t);

ics := ph(0) = 1, pc(0) = 0, pa(0) = 0, prj(0) = 0, prd(0) = 0, pgd(0) = 0, pf(0) = 0;

VARS:=map2(op,1,lhs~([sys_ode])); #The functions in the order of the derivatives on the left hand sides.
LinearAlgebra:-GenerateMatrix(rhs~([sys_ode]),VARS): A:=%[1];
# Now your system is in vector form:
diff~(Vector(VARS),t)=A.Vector(VARS);
#The eigenvalues of A are the 7 roots of p:
p:=LinearAlgebra:-CharacteristicPolynomial(A,lambda);
# Way too complicated in general.
#So we take a concrete case: All parameters are set to 1.
#But obviously you could just set them to what you want.
params:=indets({sys_ode},name) minus {t} =~1;
eval(p,params);
fsolve(%=0,lambda,complex); #The eigenvalues
#In order to find the solutions we need eigenvalues and eigenvectors:
L,V:=LinearAlgebra:-Eigenvectors(evalf(eval(A,params)));
#Notice that VARS and ics come in the same order:
VARS;
ics;
# The initial value vector:
ICS:=Vector(rhs~([ics]));
# The solution from linear algebra:
sol:=V.LinearAlgebra:-DiagonalMatrix(exp~(t*L)).V^(-1).ICS;
# We know a priori that all the solutions will be real valued functions. We would like to see that they actually are:
evalc(sol);
SOL:=fnormal~(%);
plot(SOL,t=0..5);
# Testing that our solution is actually a solution:
odetest(VARS=~SOL,eval({sys_ode},params));
fnormal(%);
simplify(%);
# If you happen to want the integral of 1-pf(t) from t to infinity you will see that it must be infinite:
int(1-SOL[-1],t=0..infinity);
##############

From your uploaded worksheet it appears that
params:={yc=1/3600,yh=1,yrd=1/180,ygd=1/180,yrj=1/180,urd=0.8,ugd=0.8,ya=1/180};

This won't change the conclusion, however. You may want to set Digits:=15 and plot from t=0..100, but otherwise the code should not need change.

You have list brackets around your expression assigned to f_jp.

It should just be

f_jp := Re(Lf((b+I*li[j]+I*evalf(Pi)*(p/m-1))/Delta));

If you want to keep the big O then you can do:

u1:=taylor(u(x-h,y),h,4);
u2:=taylor(u(x+h,y),h,4);
series(u1+u2,h);

If on the other hand you really don't want the big O anywhere, you can use mtaylor instead:

u1:=mtaylor(u(x-h,y),h,4);
u2:=mtaylor(u(x+h,y),h,4);
u1+u2;

Finally, the shortest possible versions of both:
taylor(u(x-h,y)+u(x+h,y),h,4);
mtaylor(u(x-h,y)+u(x+h,y),h,4);

@Christopher2222 Isn't the sign on U wrong?  (Note: sin instead of cos, se remark below).
If you correct that you could use events to make the ball bounce back at theta = 0:

sol := dsolve({EL, ini},numeric,events=[[theta(t)=0,diff(theta(t),t)=-diff(theta(t),t)]]);

If you add the range argument range=0..10 to dsolve and use refine=1 in odeplot, you get a nicer curve.


There is a similar example in the help page for dsolve/numeric/events.

There is no point in trying to convert the NULL response from your dsolve command to Bessel.

Edited:
Try this:
dsolve({eq1,eq2,eq3}); #General solution: involves BesselJ and BesselI
#We try to take it stepwise:
dsolve({eq1,theta(a)=A,theta(0)=1}); #Doesn't work. Probably because of BesselI(0,k*y) having a singularity at y=0).
#So try using b instead of 0 and taking limit.
res1:=dsolve({eq1,theta(a) = A,theta(b)=1});
TH:=limit(rhs(res1),b=0);
#However, notice that because of the singularity at y=0 of BesselI(0,k*y) merely imposing the condition that theta(0) be finite forces the constant _C2 in the result of
dsolve(eq1);
to be 0.
Thus although TH satisfies theta(a)=A it will not necessarily satisfy theta(0)=1.
In other words you cannot impose more than the one condition at y=a if you ask for a solution that is finite at y=0.

There are several problems:

1. A non-matching parenthesis ) is revealed when copying and pasting your dsys. I removed one of the 5 in ))))).
2. Curly braces {} are used around -1 in what should be z^(-1) in two places.
3. The initial conditions are given at z=0: Not so good because of z^(-1).
4. To solve the system you need to use a DAE method, e.g. method=rkf45_dae.
5. You can have only one condition on u.
6. The DAE-method will differentiate the second ode. The initial conditions must satisfy the second ode (before differentiation).

As an example this works:
restart;
sys := {(1-4*(diff(ln(v(z)), z)))*(diff(u(z), z))+((3/2)*z^(-1)-2*(3* z^(-1) *(diff(ln(v(z)), z))+2*(diff(ln(v(z)),z,z ))))*u(z) = 0, -z*(diff(v(z), z))-v(z)+v(z)^(1/2)*u(z) = 0};
ics:={u(0.1)=1,v(0.1)=1,D(v)(0.1)=0};
res:=dsolve(sys union ics, numeric,method=rkf45_dae);
plots:-odeplot(res,[[z,u(z)],[z,v(z)]],0.1..3);
#################

As for the analytic solution you give I find no requirements needed such as z<0 and z>-4.
sol:= {u(z)=(1+z/4)*exp(z/8), v(z)=exp(z/4)};
odetest(sol,sys);
simplify(%) assuming real;
diff~(sol,z);
convert(%,D);
eval(%,z=0);
eval(sol,z=0);

To find parameter values for which your trial solution is actually a solution, you can use solve/identity as follows.

Begin by defining odesys as you do. Change gamma in your trial solution puta to gamma1 because gamma is Euler's constant in Maple.

Then do:
res := op(odetest(puta, odesys));
SOL:=solve(identity(res=0,r),[a,alpha,beta,gamma1,lambda]);

In Maple 2015 we find 8 different solutions for the parameters.

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