@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);

@torabi I haven't answered since I have no results to show.
Let me use the opportunity to ask you a question:
How did you come up with that rather elaborate approximate solution?
I would think that the second branch of omega2 would correspond to a solution where s(x) had one zero between 0 and 1, not two.
You used the same approximate solution for the first omega2 branch, which actually corresponds to solutions with no zeros for s(x) in the open interval 0..1. (See animation at the bottom):
I tried the following revised version leading to the first branch. The point is that the approximate solution is extremely simple and leads to the same graph as the version where I used your approximate solution. To be able to see the graphs of s(x) for all the values of s1 I changed the code slightly.
Sigma:=Vector(datatype=float);
S1:=Vector(datatype=float);
AP0:=[omega2 = 0.5e-1, s(x) = 0.001*x,g(x)=0];
AP:=AP0:
i:=0:
for s1 from 0.001 to 0.06 by 0.001 do
   try
     res1[i+1]:=dsolve(dsys4 union {s(1) = s1}, numeric, initmesh = 1024, approxsoln = AP, abserr = 0.1e-4,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)]); #Same as above, no change.
plots:-display(seq(plots:-odeplot(res1[i],[[x,s(x)],[x,100*g(x)]],caption=typeset(s(1)=i*0.001)),i=1..nops([indices(res1)])),insequence);

From a very limited test I conclude that the prime notation is only used when lower case x is used. If lower case t is used then the dot notation is used.

These are not easily readable:

seq(diff(f(t),t$i),i=1..9);

This is much better:

seq(diff(f(x),x$i),i=1..7);

@vv Thanks for clearing up my confusion.

@vv I'm still confused about your use of the phrase "no collocation". Your points are also "special": They are chosen not to include x=0 and be equidistant.

@vv For the simple example I gave I found:
Collocation at Chebyshev zeros (mapped to [0,1]) with number of points = n+1 = 9, (n=8), and a basis consisting of the n+2 functions [seq(x^(i/2),i=0..n+1)] (and using fsolve) gave a maximal error of 7*10^(-9), whereas collocation with the n+1 equidistant points [seq(i/(n+1),i=1..n+1)] I got a maximal error 1.2*10^(-7).
Note: Your approach uses collocation too. Collocation just means that you ask the desired equation (in your case F(x)=0) to be satisfied (exactly) at a finite number of points. This leads to the same number of equations, which (together with boundary conditions) can be solved in various ways. I use fsolve unless it fails, only then I go to LSSolve. Alternative approaches involve integration.

@vv What I meant by collocation is exhibited in your treatment in the line
evalf(add(F(k/nx)^2,k=1..nx)):
where you evaluate F at the points seq(k/nx,k=1..nx).
The difference between my way of doing it and yours is just that I use Chebyshev zeros instead (translated to the interval [0,1]) and that I most often use fsolve to find the coefficients.