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

@one man Your comment "Unclear" is ironically unclear to me.
If you give us a concrete problem, we can have a look at it.

@Thomas Richard If I do (in Maple 2016.1)

restart;
verify(cos(u)+sin(u), sqrt(2)*cos(u-Pi/4), equal);
## I get FAIL
The only special thing I have is a maple.ini with kernelopts(floatPi=false):
Changing it to true doesn't make any difference.

@Un4 In Maple 2015 this gives a reasonable plot with 10^6 instead of undefined, which works in Maple 2016:

yield_h := proc (sigma__x, sigma__y, tau__xy, sigma__cx, sigma__cy, f__45, f__cx, f__cy, tau__u, f__tx, f__ty, alpha)
  local fr:=f__r(sigma__x, sigma__y, tau__xy, f__tx, f__ty, alpha),
        fh:=f__h(sigma__x, sigma__y, tau__xy, sigma__cx, sigma__cy, f__45, f__cx, f__cy, tau__u);
if evalf(fr) < 0 then fh else 10^6 end if end proc;
yield_r := proc (sigma__x, sigma__y, tau__xy, sigma__cx, sigma__cy, f__45, f__cx, f__cy, tau__u, f__tx, f__ty, alpha)
  local fr:=f__r(sigma__x, sigma__y, tau__xy, f__tx, f__ty, alpha),
        fh:=f__h(sigma__x, sigma__y, tau__xy, sigma__cx, sigma__cy, f__45, f__cx, f__cy, tau__u);
if evalf(fh) < 0 then fr else 10^6 end if end proc;
plots:-implicitplot3d([rcurry(yield_h,5$9),rcurry(yield_r,5$9)], -10..10, -10 .. 10, 0 .. 5, style = surfacecontour, grid=[50,50,50], axes = normal,shading=[zhue,zgrayscale);


In Maple 2016 with undefined instead of the number 10^6:

You will notice the difference if you have both versions, turn the graphs and look under the skirt :-)

@Un4 By searching MaplePrimes I found

http://www.mapleprimes.com/posts/137640-Two-Variants

where Kitonum's example

F := piecewise(x >= 0 and y >= 0 and z >= 0, 2*x+3*y+z-6, undefined):
plots:-implicitplot3d(F, x = -3 .. 3, y = -3 .. 3, z = -1 .. 7, color = red, axes = normal, scaling = constrained, style = surface, numpoints = 50000, view=[-3 .. 3, -3 .. 3, -1 .. 7]);

was described by Markiyan Hirnyk as not working in Maple 16.01, but OK in Maple 13.
The example works in Maple 2016, but not in Maple 2015.