@olivertwist The best way to find out why something is needed is to try with and without.
In this case (for me at least) I had to isolate the derivatives for shoot to work, which meant passing from ODE1 to ODE2.
If I'm right, this simply means that the shoot procedure was written under the assumption that the derivatives are isolated on the left.

@Carl Love Yes, certainly. I just tried your code in Maple 8. It worked.
The copyright is
op(3,eval(`dsolve/numeric/bvp`));
    `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`
## I have also tested the code in Maple 6 and 7. `dsolve/numeric/bvp` doesn't exist in Maple 6, but does in Maple 7, and your code works in Maple 7, but not in Maple 6.

@tomleslie Yes, I missed that one.

@Doug Meade I have been doing what you describe at the end. I use dsolve with the parameters option and fsolve.
For the system in this question your shoot was roughly fifty times faster than mine. I attribute that to my use of fsolve.
Note: I just added LSSolve as a solver. That makes my procedure go about as fast as yours at least in this example.

@vv I found the horses in the subfolder 'images', so I needed:
X:=Read(cat(kernelopts(mapledir),"/data/images/fjords.jpg")):

@Carl Love Yes, I can see that I (among other people) made a comment to that question. I still have a copy.

@John Fredsted Your method works for me too. In fact it worked when I opened the next to the last section only (by left clicking on the overlapping triangles).

@vv Thanks. My (embarrassing) mistake.

@vv I may not understand what you are saying, but both limits below exist and are equal:

eval(x-a+y-2, [x=a+rho*cos(theta), y=b+rho*sin(theta)]);
limit(%, rho=0);
# and
limit(x-a+y-2,{x=a,y=b});

The latter (described in ?limit/multi) should be quite equivalent to the former in my understanding of limit.

# What confuses me in your answer is the reference to "the previous limit". What was that previous one?

 # Note: I made an embarrassing mistake. See examples by vv below.

@tomleslie I'm guessing that once the system is written as an obviously linear system fsolve will use that fact, so that there is no need for LinearSolve.
As a rather trivial example take sys and sys2 below, where sys2 is linear but sys is not. Assuming that A[i]<>0 for all i the two are equivalent, but fsolve has more trouble with sys:
restart;
sys:={seq(1/A[i]*add(j*A[j]+i*j,j=1..i)=0,i=1..3)};
sys2:={seq(add(j*A[j]+i*j,j=1..i)=0,i=1..3)};
debug(fsolve);
fsolve(sys);
fsolve(sys2);
##### Timing comparisons:
restart;
sys:={seq(1/A[i]*add(j*A[j]+i*j,j=1..i)=0,i=1..500)}:
CodeTools:-Usage(fsolve(sys)):
memory used=1.23GiB, alloc change=0.61GiB, cpu time=15.61s, real time=14.81s, gc time=11.19s
restart;
sys2:={seq(add(j*A[j]+i*j,j=1..i)=0,i=1..500)}:
CodeTools:-Usage(fsolve(sys2)):
memory used=256.07MiB, alloc change=26.71MiB, cpu time=2.39s, real time=2.41s, gc time=171.88ms
##############
To Jens I should add that when fsolve turns unevaluated (as it did for you) this doesn't mean that it hasn't performed any work. It means that it wasn't able to find a solution.

@torabi I can only report what I find. I shall leave it to you to compare the results to the results reported in the paper where you found the graphs.

@torabi The approxsoln  g(x)=0  and  s(x) = 0.001*x used in trying to find the first branch uses an s(x) with no zero in the interval (0,1). g(x)=0 is used because g will be small. But the proof lies in the pudding (as I think the saying goes): It works.
In correcting what you want plotted I'm left confused. You write:
plot(Sigma/(L*sqrt(m/E_a)), S1*L, labels = [sigma, s(1)]);

but my Sigma is the vector with components sqrt(omega2(1))/(L*sqrt(m/E_a)).
Probably (?) you meant plot(Sigma,S1*L, labels = [sigma, s(1)]); #The label s(1) should also be changed then.
Like in
plot(Sigma, S1*L, labels = [sigma, 'L'*s(1)]);

@Axel Vogt RES is talking about having a Maple reader, not Maple itself. I suppose he has the free Maple Player, but from that you cannot export as text as far as I see.

I don't know about the reader, but the task should be quite simple from Maple itself.
I don't suppose you would want to upload one (or all) of the .mws files to MaplePrimes?

@torabi I get what I would call a second branch here (notice that the loop is going backwards).

res1:='res1':
Sigma := Vector(datatype = float);
S1 := Vector(datatype = float);
AP:=[omega2 = .033, s(x) = -x^2*(0.64-x^2)/30, g(x) = -x/100]; #One zero in (0,1).
i := 0;
for s1 from 0.01 by -0.001 to 0.001 do
  try
    res1[i+1] := dsolve(dsys4 union {s(1) = s1}, numeric, initmesh = 4096, approxsoln = AP, abserr = 0.1e-3, maxmesh = 8192);
    i := i+1;
    AP := res1[i];
    Sigma(i):=eval(sqrt(omega2(1))/(L*sqrt(m/E_a)), res1[i](1));
    S1(i) :=s1
  catch:
    print(s1);
    print(lasterror)
  end try
end do;
plot(Sigma, S1, labels = [sigma, s(1)]);

plots:-display(seq(plots:-odeplot(res1[i],[[x,s(x)],[x,100*g(x)]],caption=typeset(s(1)=0.01-(i-1)*0.001)),i=1..nops([indices(res1)])),insequence);