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

@Christopher2222 You cannot have an event written like you do:

[diff(theta(t),t)=0,a:=[op(a),t]]

The assignment is clearly out of place, but a=[op(a),t] would also be wrong.
Before I try to elaborate I should like to know what 'a' is supposed to be. What has it got to do with EL. Is it supposed to be a discrete variable?

My guess is that you intended something like

sol3 := dsolve({EL, ini} union {b(0)=[]},numeric,discrete_variables=[b(t)::list],events=[[theta(t)=0,diff(theta(t),t)=-diff(theta(t),t)],[diff(theta(t),t)=0,b(t)=[op(b(t)),theta(t)]]],range=0..10);

But you get an error telling you that discrete variables must be of type boolean, integer, or float.

This restriction is explicitly mentioned in the help page
?dsolve,numeric,events
under the heading Discrete variables:

Only the types boolean, integer and float are currently available.

################
If my guess is right you may try something like this instead:

sol3 := dsolve({EL, ini} union {b(0)=0,TT(0)=0},numeric,discrete_variables=[b(t)::float,TT(t)::float],events=[[theta(t)=0,diff(theta(t),t)=-diff(theta(t),t)],[diff(theta(t),t)=0,[TT(t)=t,b(t)=theta(t)]]],output=listprocedure);
B,TTp:=op(subs(sol3,[b(t),TT(t)]));
plots:-odeplot(sol3,[[t,theta(t)],[t,b(t)],[t,TT(t)]],0..10);
L:=[]:
for tt from 0 to 10 by .5 do L:=[op(L),[TTp(tt),B(tt)]] end do:
L;
ListTools:-MakeUnique(L);

@Christopher2222 How about
solL := dsolve({EL, ini},numeric,events=[[theta(t)=0,diff(theta(t),t)=-diff(theta(t),t)]],range=0..10,output=listprocedure);
TH,TH1:=op(subs(solL,[theta(t),diff(theta(t),t)]));
fsolve(TH1,[2]);
TH(%);
fsolve(TH1,[4]);
TH(%);
## or
sol2 := dsolve({EL, ini} union {b(0)=0},numeric,discrete_variables=[b(t)::float],events=[[theta(t)=0,diff(theta(t),t)=-diff(theta(t),t)],[[diff(theta(t),t)=0,theta(t)>0],b(t)=theta(t)]],range=0..10);
plots:-odeplot(sol2,[[t,theta(t)],[t,b(t)]],0..10);
sol2(5);

@Athar Shahabinejad I see no way of doing that.
The characteristic polynomial is a polynomial of degree 6 with rather general coefficients. In fact it appears that they are so general that any monic polynomial of degree 6 can be written in that form. Thus it is well known that no formula exists (or can be found) for the 6 roots of that polynomial.
Unless you give values to all or some of the parameters I cannot see that you can get anywhere.