@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.

@Kitonum Certainly this is much nicer. I don't know what the original question was hinting at, but if it had to do with exhibiting a weakness in verify then it succeeded.

@Harry Garst In the help page for verify/equal it says:

The verify(expr1, expr2, equal) function returns true if
                signum(0, expr1 - expr2, 0) = 0
 and false if the call to signum returns a numeric object. Otherwise, it returns FAIL.


If you try

signum(0,abs(x)- sqrt(x^2), 0);

in Maple 2015 you get 
i.e. you don't get 0 nor a numeric object. Thus verify must return FAIL according to its help page.
However, if you use simplify on the signum-result then you get 0.
###################
I haven't been using verify (ever?) and cannot say that I'm impressed:
verify(I*sin(7),cos(8),equal);
signum(0,I*sin(7)-cos(8), 0);
verify(I*sin(7.),cos(8.),equal);
signum(0,I*sin(7.)-cos(8.), 0);

@Athar Shahabinejad I think I understood what you meant. The problem is that you cannot find a formula for the eigenvalues and thus not for the eigenvectors either. Therefore you cannot find a formula for the solution to your odes. Hence no formula for pf(t) and consequently you cannot evaluate the integral.