@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

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