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

@vv The following ide has a simple solution which is not C^1 on the closed interval [0,1].
ide := x^2*diff(y(x), x, x)+2*x*diff(y(x), x)-3/4*y(x)-int(y(t), t = 0 .. 1);
##
A solution is y(x)=x^(1/2)-8/21 and all constant multiples of that.
Using y(0)=-8/21 and ide with the same type basis as I used above [seq(x^(i/2),i=0..n+1)] and Chebyshev collocation points I get very good agreement between numerics and the exact solution.

@Kitonum I was able to reproduce your second solution by the same method as used by vv, but with a different basis.
I found
 y(x) = 1.+.44556351*sqrt(x)-10.06038868*x+56.55938008*x^(3/2)-192.5721398*x^2+429.3716449*x^(5/2)-631.4150415*x^3+610.7682516*x^(7/2)-373.7890882*x^4+131.3990087*x^(9/2)-20.2071904*x^5

by using the basis [seq(x^(i/2),i=0..n+2)] (n=8) with collocation at the n+1 Chebyshev zeros given by
[0.75961235e-2, 0.669872980e-1, .1786061952, .3289899286, .5000000000, .6710100714, .8213938048, .9330127020, .9924038765].
Of course I included y(0)=1 and y(1)=1.5 as additional equations.

@vv I have no intention of writing a paper on this, but do try to treat the problem seriously. Otherwise this is all a waste of time.

I find it interesting (disturbing?) that the solution found by Kitonum's approach with y(0)=1, y(1)=1.5 actually seems to persist as n is increased, while attempts at finding a solution by the method you described (one that I have also been using for a while) doesn't seem capable of producing a similar result.
Note: See correction below.

Why is this an emergency?

@vv Your argument is correct if y2'(x) and y2''(x) are bounded in an interval (0,eps), eps>0. That may not be the case, however. Just consider the ode
odeH:=x^2*diff(y(x), x, x)+50*x*diff(y(x), x)-35*y(x) = 0;
## It has the general solution
y(x) = _C1*x^(-49/2+(11/2)*sqrt(21))+_C2*x^(-49/2-(11/2)*sqrt(21)).
Only the trvial solution y(x)=0 satisfies the boundedness requirement above since one exponent is negative (~-49.7) and the other is positive, but less than 1 (~ 0.7).
This doesn't prove that y2 doesn't satisfy the boundedness requirement, but it suggests that we should be careful here.
It could be that x^2*y2''(x)+50*x*y2'(x) had a finite nonzero limit as x -> 0 (from the right).