@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

@vv You are replacing y(t) in the integral by p defined by p:=add( a[k]*x^k,k=0..n), thus you might just as well begin by replacing y(t) by y(x) in the integral after which the integration can be performed without knowing y. Thus we have an ordinary ode with a regular singular point at 0.

ode:=x^2*(diff(y(x), x, x))+50*x*(diff(y(x), x))-35*y(x) - (1-exp(x+1))/(x+1)-(x^2+50*x-35)*exp(x)-int(exp(x*t)*y(x),t=0..1);
ode0:=select(has,ode,y(x));
DEtools:-indicialeq(ode0,x,0,y(x));
solve(%,x);

#But why is it reasonable to replace y(t) by y(x)?

To answer my own question I tried the following, where I just repeat your lines to begin with.

restart;
ide:=x^2*(diff(y(x), x, x))+50*x*(diff(y(x), x))-35*y(x) - (1-exp(x+1))/(x+1)-(x^2+50*x-35)*exp(x)-int(exp(x*t)*y(t),t=0..1):
n:=8: nx:=8:
p:=add( a[k]*x^k,k=0..n):
F:=unapply( eval(ide, [y(x)=p,y(t)=p]),x);
d:=add( F(k/nx)^2,k=1..nx );
s:=Optimization:-NLPSolve(d,iterationlimit=10000);
yy:=eval(p,s[2]);  # approx solution
#numapprox[infnorm](eval(F(x),s[2]),x=0..1); #This just checks the approximate F(x) not ide.
##The following I do in two steps to avoid loosing the kernel. I shall report that as an SCR.
Y:=unapply(yy,x);
#eval(ide,y=unapply(yy,x)); #Looses the kernel
numapprox[infnorm](eval(ide,y=Y),x=0.1..1); #Using interval x=0.1..1 (0..1 is unfair).
plot(eval(ide,y=Y),x=0..1,-2..2);

@Kitonum As pointed out by vv the behavior is documented:
"In the case of infinite sums, Levin's u-transform is used, which has the additional effect that sums that formally diverge may return a result that can be interpreted as evaluation of the analytic extension of the series for the sum (see the examples below)."

You could try:
restart;
gg := y-> Sum(2^jj*y^jj, jj = 0 .. infinity);
debug(`evalf/Sum/LevinU_Symb/1_inf`);
evalf(gg(2));

@torabi You can do as follows.
The idea is to start with the approximate solution provided by you. But in the next step the previously found solution is used for the new value of s1. For convenience the results are kept in vectors.

Sigma:=Vector(datatype=float);
S1:=Vector(datatype=float);
AP:=[omega2 = 0.5e-1, s(x) = (56330/77299)*x^6-(25403/9094)*x^5+(1214873/309196)*x^4-(182218/77299)*x^3+(1/2)*x^2, g(x) = -(32119348934425/2772462826064)*x^5+(133402088597045/5544925652128)*x^4-(35403636020221/2772462826064)*x^3-(1/2)*x^2+x];
i:=0:
for s1 from 0.001 to 0.06 by 0.001 do
   try
     res1:=dsolve(dsys4 union {s(1) = s1}, numeric, initmesh = 1024, approxsoln = AP, abserr = 0.1e-4,maxmesh=8192);
     AP:=res1;
     i:=i+1;
     Sigma(i):=eval(sqrt(omega2(1))/(L*sqrt(m/E_a)),res1(1));
     S1(i):=s1;
   catch:
     print(s1);
     print(lasterror);
   end try;
end do;  
plot(Sigma,S1,labels=[sigma,s(1)]);

@Abdoulaye See my answer in
http://mapleprimes.com/questions/213864-Help-For-Plotting-

@Les In a Windows system a place to put the file (name maple.ini) is
C:/Users/xxxx/maple.ini
where xxxx should be replaced by your user name on the computer.
The file may not exist, but just use any text editor to create it.

@torabi The reason is that I used ZZ instead of what you actually had: Z^2.
So change pxx, pyy, pzz like this:
pxx:=subs(X=X^2,x=(x,x),eval(px));
and similarly for pyy and pzz.

But the noncommutativity issue is important and remains to be dealt with.

@torabi I made the following naive attempt not thinking of the noncommutativity of differential operators and multiplication operators. See note at bottom. I don't believe it was a problem in your previous example.

I removed the global declaration in px entirely. It is not necessary anyway.

indets(B,name); #We find (besides x,y,z, and pi) X, XY, XZ, Y, Z, ZY.

px := proc (u1) local u;
  u := expand(u1);
  if type(u, `+`) then return map(procname, u) end if;
  if not type(u, `*`) then return u end if;
  if member(X, {op(u)}) then diff(u/X, x) else u end if
end proc;
py := subs(X = Y, x = y, eval(px));
pz := subs(X = Z, x = z, eval(px));
pxx:=subs(X=XX,x=(x,x),eval(px));
pyy:=subs(X=YY,x=(y,y),eval(px));
pzz:=subs(X=ZZ,x=(z,z),eval(px));
pxy:=subs(X=XY,x=(x,y),eval(px));
pxz:=subs(X=XZ,x=(x,z),eval(px));
pzy:=subs(X=ZY,x=(y,z),eval(px));
P:=`@`(px,py,pz,pxx,pyy,pzz,pxy,pxz,pzy);
R:=map(P,B);
###############
### Time for a very critical warning: Differential operators don't commute with multiplication operators.
In your recent worksheet Q has e.g. this element:
Q[2,1];
## In your input you write it as
2/(Pi*a*(a*y+b))*ZY;
## Ordinary multiplication in Maple is commutative, so that could be a problem depending on what your actual meaning was.
This operator will be handled by pzy as
ZY*2/(Pi*a*(a*y+b)); #In that order and was that intended?
## Now compare
diff(2/(Pi*a*(a*y+b))*f(x,y,z),z,y);
## with
2/(Pi*a*(a*y+b))*diff(f(x,y,z),z,y);
## My conclusion is that you must rewrite Q so that the noncommutativity is respected. Maybe use &* or . (matrix multiplication).
The whole must be thought over and another approach used. The problem is not that Q and N are matrices.
You must first be able to handle an operator like 2/(Pi*a*(a*y+b))&*ZY, where ZY stands for diff(..., y,z). I only used &* here to keep the order if entered in Maple.

@mathematicianS In your original formulation of the boundary conditions you have f'''(0), yet the ode for f is only of order 3. You must reformulate that so that only orders lower than the highest derivative order in the ode appears. 
In your first order systems version the same problem appears since D(v)(0) appears, i.e. v'(0).

Closely related questions have been asked and answered often in the past. For anybody to help you give us the odes and the boundary values as text so we don't have to type.
Since solution will be numerical we shall need the values of lambda, Pr, s, sigma, a, and b.

@vv In a programmatic context it appears that the repeated creation of the same procedure keeps the address.
In these examples I removed option remember since it seems to be irrelevant for my point:

restart;
F:=proc() 0 end proc;
addressof(eval(F));
F:=proc() 0 end proc;
addressof(eval(F)); #Different from above
to 10 do F:=proc() 0 end proc; print(addressof(eval(F))) end do: #All addresses are the same
to 10 do F:=proc() 0 end proc; print(addressof(eval(F))) end do: #All the same, but different from just above
F:=proc() 0 end proc; addressof(eval(F)); #Try executing this very line twice. Result: Different addresses.
p:=proc() global f; f:=proc() 0 end proc; NULL end proc; #global might as well have been local, no difference
p();
addressof(eval(f));
p();
addressof(eval(f)); #Same as just above
############
These examples remind me of an unrelated discussion of a weird bug that only occurred when the user entered the code. I have forgotten any details.

@Carl Love Rather strange.
I made this observation. Compare the following 2 procedures.
q1 has your f, q2 has a modified version of your f. f in q2 just has an evaluation of a variable k local to q2.
q1 apparently needs forget, q2 doesn't.
Note added: If the NULL output is changed to f, then it appears that the two successive outputs from q1() evaluate to the same procedure, whereas this is not the case for q2.
Conclusion: The behavior of q1 is not a bug. The procedure with the local name f (name new at its invocation) is the same in each call. The memory table belongs to the procedure, not to its (many) names. Thus the memory table remains.
But the procedure f produced by q2 at each call is so to speak "parametrized" by the local variable k. Thus these procedures f are all different.
####
restart;
q1:=proc() local f, k; #k not used anywhere though
   f:=proc() option remember; 0 end proc;
   print("Entering:", op(4,eval(f)));
   f(89);
   print("Leaving:",op(4,eval(f)));
   NULL
end proc;
q1();
q1();
q2:=proc() local f,k;
   f:=proc() option remember; k; 0 end proc;
   print("Entering:", op(4,eval(f)));
   f(89);
   print("Leaving:",op(4,eval(f)));
   NULL
end proc;
q2();
q2();