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.

Experimenting just a little with changing (i,j)->`if`(j>=i,1,0) to e.g.
(i,j)->`if`(j>=i,j,0) seems to show that the first nonzero element is returned by add in the case of testM.

Trying seq instead of add seems to work as expected. For testM the nonzero elements are taken first.

Trying mul instead of add gives 0 as expected for nm, but the same as add for testM, i.e. the first nonzero element.

@Carl Love Yes, add(x, x= myMat[5, 1..-1]) works as expected in Maple 16. However, it doesn't in Maple 17, 18, and 2015.

@Carl Love Your version add(x, x= myMat[5, 1..-1]) doesn't work (for me).
However, add(myMat[5,i],i=1..40) does.