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

@matteoziff To me the most natural thing to do would be to rewrite the ode as a first order system (sys below) in the standard way using the new additional dependent variable y given by y(t) = diff(x(t),t). Then after that use DEplot from the DEtools package to make a direction field

DEtools[DEplot](sys,[x(t),y(t)], t=0..1, x = -2.1..2.1, y = -1500..1500);

and show that together with the plot from odeplot using display.

But it seems that you don't want to do that?

@Noreen cute I tried the worksheet in Maple 12 (I don't have a working version of Maple 13).

The problem is with fsolve. The syntax used works in Maple 16, but not in Maple 12. I changed that.

So here is a revised version, which works in Maple 12 and I'm sure also in Maple 16.

MaplePrimes12-08-23.mw

@Noreen cute I tried the worksheet in Maple 12 (I don't have a working version of Maple 13).

The problem is with fsolve. The syntax used works in Maple 16, but not in Maple 12. I changed that.

So here is a revised version, which works in Maple 12 and I'm sure also in Maple 16.

MaplePrimes12-08-23.mw

Since this problem is quite similar to the one you posed yesterday:

http://www.mapleprimes.com/questions/136745-How-To-Calculate-Dual-Solution-Problem

I tried the same method. It works fine with b = 6, but certainly failed for b = 20.

But the direct approach you are using also works with b = 6 (even b=6.5). When Pr = 2 you may need something like
sol:=dsolve(eval({de1},Pr=2),numeric,continuation=k);

which worked for me. It starts solving the problem with k=0 and then works its way up to k=1.

@Noreen cute What kind of error?

You could try to execute the following worksheet in Maple 16:

MaplePrimes12-08-23S.mw