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

@Kitonum To your question: I don't know. I haven't been able to reproduce the second solution by an approach similar to the one used by vv, but including boundary conditions y(0)=1, y(1)=1.5.
My attempts at proving uniqueness of the problem with y(0) = 1 also has failed so far beginning like this:

Your second solution y2 (if it is a solution) satisfies y2(0)=1 and y2(1)=1.5. Furthermore it appears that it is less than y1(x) = exp(x) on the half-open interval (0,1]. Let us assume that then.
Thus y(x) = y1(x) - y2(x) > 0 on the interval (0,1].

The difference y satisfies the equation ("=0" assumed as always):
ideH:=select(has,ide,y); # where I use the notation ide used by vv.
## Now let us replace the integral int(exp(x*t)*y(t),t=0..1) by phi(x):
ode:=subsop(-1=-phi(x),ideH);
odeH:=eval(ode,phi=0);
## Considering phi known the corresponding homogeneous equation odeH has the general solution:
dsolve(odeH);
        y(x) = _C1*x^(-49/2+(11/2)*sqrt(21))+_C2*x^(-49/2-(11/2)*sqrt(21))
## One of the exponents is negative, which is what I hoped to be able to use to prove uniqueness.
##ode can be solved for y in terms of phi by dsolve (not using y(0)=0 since the outcome of that is problematic).
## But I'm stuck.

@John Fredsted Yes, I agree.

@John Fredsted While your objection is called for here, it is too strong.

Example:
restart;
p:=proc(x) x:=8 end proc;
p(a);
a;
p(5); #Error as the number 5 cannot be assigned to
p(a); #Error as the variable a evaluates to the number 8 before p acts on it.
# To avoid the latter error we can do:
q:=proc(x::uneval) x:=67 end proc;
q(c);
c;
q(c);

@ThU Yes post the code.

The following would produce the error (but won't mention assuming):

restart;
t0:=7,8; #Notice the comma instead of the dot it ought to be.
F:=1: w:=2:
evalf(F*cos(w*t0));
   
Error, (in cos) expecting 1 argument, got 2