@marram I have edited my original answer quite a lot since it was not correct. I looked up the paper you referred to and now know what the authors call the characteristic equation in this context. Sorry about the confusion!

Be aware that in the Maple version they used in 1997 the LinearAlgebra package was not yet introduced. It came with Maple 6. The old package linalg should not be used these days.

@marram I have edited my original answer quite a lot since it was not correct. I looked up the paper you referred to and now know what the authors call the characteristic equation in this context. Sorry about the confusion!

Be aware that in the Maple version they used in 1997 the LinearAlgebra package was not yet introduced. It came with Maple 6. The old package linalg should not be used these days.

@J4James Since Nc appears in the boundary condition at 0 you need to include that in the eval part, like this:

dsolve(eval({eq2,theta(9)=0,D(theta)(0)=-Nc*(1-theta(0))},{beta=1,Ec=0.5,Q=-0.5,Nc=5,pr=3,f(eta)=F(eta),diff(f(eta),eta,eta)=F2(eta)}),numeric,known=[F,F2]):

@J4James Since Nc appears in the boundary condition at 0 you need to include that in the eval part, like this:

dsolve(eval({eq2,theta(9)=0,D(theta)(0)=-Nc*(1-theta(0))},{beta=1,Ec=0.5,Q=-0.5,Nc=5,pr=3,f(eta)=F(eta),diff(f(eta),eta,eta)=F2(eta)}),numeric,known=[F,F2]):

It would be helpful if you gave us the system.

@J4James This is not simple. I tried and have uploaded a worksheet which does get something looking like fig. 2a in the paper. Since infinity = 9 it is not surprising that there is a substantial difference. However, the general shape is good. As could be expected it is tricky at the leftmost end.

MaplePrimes12-09-06B.mw

Although you don't give me all the necessary details, the loop is not OK:

A := []; for e from -1 by .1 to 3 do solk := dsolve({icsys, op(subs(k = 10^e, sysk))}, numeric,

 known = [A, A1, A2, A3]); deltk := subs(solk, Delta[d](t)); 

A := [op(A), [e, log10(deltk(.9)^2)]] end do;

A is defined to be a list initially, thus not the name of a procedure. That alone would cause an error.
Secondly A is again defined to be a list at the end of the loop and furthermore

[e, log10(deltk(.9)^2)] is again a list, and not even of names of procedures.

Although you don't give me all the necessary details, the loop is not OK:

A := []; for e from -1 by .1 to 3 do solk := dsolve({icsys, op(subs(k = 10^e, sysk))}, numeric,

 known = [A, A1, A2, A3]); deltk := subs(solk, Delta[d](t)); 

A := [op(A), [e, log10(deltk(.9)^2)]] end do;

A is defined to be a list initially, thus not the name of a procedure. That alone would cause an error.
Secondly A is again defined to be a list at the end of the loop and furthermore

[e, log10(deltk(.9)^2)] is again a list, and not even of names of procedures.

@J4James Is Lambda (or Omega) given in advance and the parameters s and lambda to be determined? Or is s also given and only lambda to be determined?

Do you really mean H(0) versus lambda1? Isn't H(0) = 1 during the whole process?

restart;
p:=plots:-implicitplot(sqrt(x^2+y^2)-sqrt((x-4)^2+(y-3)^2) = 5, x = -20 .. 20, y = -20 .. 20, numpoints=10^6,
 thickness = 5, scaling = constrained):
dt:=plottools:-getdata(p):
M:=<seq(dt[k,3],k=1..nops([dt]))>:
plot(M,thickness=5);

You can add the argument style=point. Then at least the points found are shown.

The disappearance of the the two lines satisfying the equation for values less than 5 can be seen in this animation:

animate(implicitplot,[sqrt(x^2+y^2)-sqrt((x-4)^2+(y-3)^2)=a, x =-20 .. 20, y = -20 .. 20, gridrefine=2],a=4.9..5);

Is it not rather a weakness than a bug?

Before answering this I tried the Subscribe link to your question. It seemed to work. At this moment it says "Subscription Created".

@Amare dsolve/numeric doesn't return a module.
If sol2:=dsolve(....,numeric); then by default a procedure is returned (a procedurelist).
If sol2:=dsolve(....,numeric,output=listprocedure); then a list of procedures is returned.

@Amare dsolve/numeric doesn't return a module.
If sol2:=dsolve(....,numeric); then by default a procedure is returned (a procedurelist).
If sol2:=dsolve(....,numeric,output=listprocedure); then a list of procedures is returned.