@Alejandro Jakubi Thanks for digging deeper into this. The difference between pp and qq is illuminating.

I submitted a bug report about convertsys giving the example diff(x(t),t,t) + diff(a(t),t,t) = 0.

Although surely the procedure was intended primarily for dsolve/numeric, it is not mentioned in the documentation that known functions can contain no derivatives.
I ran into this problem several years ago as mentioned above, but cannot remember if I filed a bug report at that time.

@StefKi Your system is linear in the second derivatives, so it ought to be simple to solve for those (when possible).
convertsys has to do that.

I tried doing it myself:

res:=solve(Sys,diff(GenKo,t,t)):
nops(res); #7 as expected
length(res); #Huge so do indeed use colons
FirstOrdSys := convertsys(res, [], GenKo, t, Y, YP ): #working on res instead of Sys
whattype(FirstOrdSys); #list
nops(FirstOrdSys); #4 as expected
length(FirstOrdSys); #Huge
length(FirstOrdSys[1]); #The bulk of the length
FirstOrdSys[2]; #The 'translation'
FirstOrdSys[3..4];

The total time used was less than 2 minutes.

Finally, your system doesn't contain derivatives of known functions (i.e. functions that are not dependent variables). Therefore the handling (or lack thereof) of such systems is irrelevant for your system.

@momoklala I have no problem with S:=-0.5.  I get

[-1.18614066271840, -1.17788387531395, -1.16926687825642, -1.16025994699062, -1.15082789257908, -1.13052356717881]

I use Maple 2015. I shall test it in Maple 18 later.
"Later": Well, it gives me the exact same results without any problem.

So what are you doing?

Several years ago I wrote a KonvSys procedure tailored for teaching purposes. I simply built it on top of `DEtools/convertsys`. Looking at it now I see that I ran into the very same problem raised here.
I solved it at that time by using the following modification, which still works:

STR:=convert(eval(`DEtools/convertsys`),string);
STR1:=convert('indets(ndeqns,'specfunc(anything,diff)')',string);
STR2:=convert('select(has,indets(ndeqns,'specfunc(anything,diff)'),u)',string);
konvsys:=parse(StringTools:-Substitute(STR,STR1,STR2));

Now just use konvsys instead of `DEtools/convertsys`.

@Alejandro Jakubi I believe you are wrong. Changing the local to something else like aa doesn't change the behavior. I maintain that the difference in behavior between using 'a' or 'p' is due to the ordering I mentioned.

Consider the following extremely stripped version of the procedure `DEtools/convertsys`. The problems appear early so it is not very long:

restart;
pc:=proc(deqns,nu)
local diffset,a,diffmax,idiff,j,ndeqns;
if type(deqns,`=`) then ndeqns:={deqns} elif type(deqns,'list') then ndeqns:={op(deqns)} else ndeqns:=deqns end if;
diffset:=indets(ndeqns,'specfunc(anything,diff)');
for a in nu do diffmax[a]:=0 end do;
printf("%a\n",Diffmax[a]=diffmax[a]);
printf("%a\n",convert(diffset,list));
for idiff in diffset do
  j:=0;
  printf("%a %a\n",'j'=j,idiff1 =idiff);
  while type(idiff,'specfunc(anything,diff)') and nops(idiff)=2 do
    j:=j+1;
    idiff:=op(1,idiff);
    printf("%a %a\n",'j' =j, idiffwhile = idiff)
  end do;
  printf("%a %a\n",idiff2=idiff,'op'(0,idiff)=op(0,idiff));
  diffmax[op(0,idiff)]:=max(j,diffmax[op(0,idiff)]);
  printf("%a\n",Diffmax[op(0,idiff)]=diffmax[op(0,idiff)])
end do;  
end proc:
##The two test equations:
test:=diff(x(t),t,t)+diff(a(t),t,t)=0;
testp:=diff(x(t),t,t)+diff(p(t),t,t)=0;
#Testing
pc(test,[x]);
pc(testp,[x]);
pc(test,[a]);
pc(testp,[p]);
##
## I'm convinced from this that the problem is as simple as exhibited right here in procedure pp:
pp:=proc(i,j) local dm; dm[a]:=max(i,dm[a]); dm[a]:=max(j,dm[a]) end proc;
pp(1,1);
pp(1,2);

#The problem surfaces in  pc(test,[x]) and not in the other 3 tests because of evaluation rules in procedures.

Notice that if you wrap max(j,diffmax[a]) in eval, then all 4 tests produce the error.

@StefKi In order to answer this we would have to see your system. You could upload a worksheet.

@StefKi What I meant by calling a(t) a known function as opposed to the unknown x(t) was that in your call to convertsys you are giving x(t) as the unknown function.

If a(t) is actually an unknown function and if you have several odes then just use convertsys on the system including a(t) as an unknown as in this simple example:

ode1:=diff(x(t),t,t)+diff(a(t),t,t)=0;
ode2:=diff(a(t),t,t)+x(t)*diff(a(t),t)=g(t);
DEtools[convertsys]({ode1,ode2},[],[x(t),a(t)],t,Y,YP);
#################
If your situation is not as above i.e. a(t) is not an unknown in that sense, then freezing and thawing might be what you need.
Again a simple example:

restart;
ode1:=diff(x(t),t,t)+diff(a(t),t,t)=0;
ode2:=diff(y(t),t,t)+x(t)=diff(b(t),t)+g(t);
vars:=[x(t),y(t)]; #The unknowns
AB:=remove(has,indets({ode1,ode2},specfunc(anything,diff)),vars);
ABeqs:=AB=~freeze~(AB);
thaw(ABeqs); #Just testing the thawing.
DEtools[convertsys](subs(ABeqs,{ode1,ode2}),[],vars,t,Y,YP);
thaw(%);
############
You might ask why some device like this freezing-thawing process is not part of the procedure convertsys.
My guess is that convertsys was primarily intended for numerical solution of odes, where conversion to a first order system is routinely done and where you will not encounter 'known' functions that are actually not known as is the case with your a(t).

@momoklala Try this:

restart;
Digits:=15:
Eq1 := diff(f(eta), eta, eta, eta)+f(eta)*diff(f(eta), eta, eta)-diff(f(eta), eta)^2 = 0;
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;
q:=proc(n) local res;
  res:=dsolve(eval({Eq1, Eq2, bcs1, bcs2}, {Pr = 7,N=n}), numeric, output = listprocedure,initmesh=2048);
  subs(res,theta(eta))(0)
end proc;
plot(q,0.1..1.89,labels=[N,theta(0)]);
plot(q,0.1..2,adaptive=false,numpoints=20,labels=[N,theta(0)]);

Why do you say those values are wrong? What is your reason?

Incidentally, there is a missing assignment symbol (:=) in the assignment to bcs2.
Also it is not a very good idea to set Digits to 5. The default 10 should be fine, but if you would like to change then you should rather raise Digits to 15. However, that is not going to get you the values you (for some reason) want.

@momoklala What I meant was very simple.

Remove the line N:=1. This line assigns the value 1 to N. You obviously don't want that since you want to vary N.
Wthether you assign to Pr or nor is irrelevant, but I chose not to. I used
eval({Eq1, Eq2, bcs1, bcs2}, {Pr = 7,N=n}),
where Pr is replaced by 7 before dsolve is set to work.

@momoklala You can mimic the idea from my solution procedure p above.

Changes: Do not assign to N. No need to assign to Pr.
Here I use Pr=7 inside the procedure q.

q:=proc(n) local res;
  res:=dsolve(eval({Eq1, Eq2, bcs1, bcs2}, {Pr = 7,N=n}), numeric, output = listprocedure,initmesh=2048);
  subs(res,theta(eta))(0)
end proc;
plot(q,1..1.8); #It takes a while so you might want this one first:

plot(q,1..2,adaptive=false,numpoints=20);

Trying to evaluate q at 1.9 or higher shows that there are problems:
q(1.9);

@momoklala I have introduced a procedure p which simply computes theta(0) as a function of Pr.

It does not use your 5 values, but uses any needed by plot.

It is not at all clear to me what you want to plot. For each value of Pr you have a solution for theta(eta) for eta in 0..blt.

If you just want to plot the 5 graphs of theta in one coordinate system then it is easy:
plot([Y11,Y12,Y13,Y14,Y15],0..blt);

but that may not be what you mean.

Incidentally, in your code you have with*plots. You probably meant with(plots);

In your example what are the x and y limits?
You write "with bounds a= -sqrt(a^2-y^2) b= 0 c= 0, d= 1".
The first "bound" in particular is confusing.

@tomleslie To answer your question: I don't think anybody finds this acceptable. It is an obvious bug.