There is no polynomial solution to your integral equation. The solution is u(x)=sinh(x).
But I guess you know that. Thus you must cut off. In your case leave out the x^4 term:

restart;
intE:=[diff(u(x), x) = 1+int(u(t), t = 0 .. x), u(0) = 0]; #I assume this was your intE (basically)
ode:=diff(intE[1],x);
dsolve({ode,u(0)=0,D(u)(0)=1}); #Answer u(x)=sinh(x) (in terms of exp)
##
U3:=unapply(add(a[i]*t^i,i=0..3),t);
intE1 := eval(intE,u=U3);
eval((lhs-rhs)(intE1[1]),intE1[2]);
pol:=collect(%,x);
seq(coeff(pol,x,i),i=0..4);
solve({%},{seq(a[i],i=1..3)}); #No solution
## Now cut off
seq(coeff(pol,x,i),i=0..3); #Forgetting  the x^4 term
solve({%},{seq(a[i],i=1..3)});
##Using dsolve:
dsolve({ode,u(0)=0,D(u)(0)=1},u(x),series,order=4);
convert(%,polynom);

Here is an example:

restart;
M:=Matrix([[f(t)*g(t),f(t)/g(t)],[f(t),g(t)]]);
inits:={f(0)=5,g(0)=2,D(g)(0)=3,D(f)(0)=-2};
diff~(M,t); #Elementwise ~ needed here
eval(%,t=0);
convert(%,D);
eval(%,inits);

I reformulated your code some.
Revised:
I use D1 instead of D since D is used for differentiating a function as in D(sin), which is cos.
Note: You forgot to multiply by k, k^2, k^3 in eq2,eq3,eq4. Since k=0 is somewhat problematic anyway, so if k<>0 it can just be forgotten from those equations.
restart;
Digits:=15:
b := 0.17e-1;
a := 0.8e-1;
n := 1;
R := r->A*BesselJ(n,k*r) + B*BesselY(n, k*r) + C*BesselI(n, k*r)  + D1*BesselK(n, k*r) ;
##Your equations:
eq1 := R(b) = 0;
eq2:=D(R)(b)=0; #First derivative of R at b
eq3:=(D@@2)(R)(a)=0; #Second derivative of R at a
eq4:=(D@@3)(R)(a)=0; #Third derivative of R at a
#The system is linear in A,B,C,D1:
M,z:=LinearAlgebra:-GenerateMatrix([eq1, eq2, eq3, eq4],[A,B,C,D1]);
#The system is homogeneous so the determinant of M must be zero if we are to have nontrivial solutions for A,B,C,D1:
det:=LinearAlgebra:-Determinant(M);
plot(det,k=0..200);
plot(det,k=0..200,view=-1..1);
#The first nonzero value for k:
k1:=fsolve(det,k=25); #Using the plot for the guess k=25
M1:=eval(M,k=k1);
LinearAlgebra:-Rank(M1); #3 as it ought to be
G:=LinearAlgebra:-GaussianElimination(M1);
G:=fnormal~(G); #Removing roundoff
LinearAlgebra:-LinearSolve(G,z,free=t);
res:=eval(%,t[1]=1);
<A,B,C,D1> = res;
Equate(<A,B,C,D1>,res);
sol1:=eval(R(r),{%[],k=k1});
plot(sol1,r=b..a);
## You can repeat the lines from the fsolve statement and down to find the next eigenvalue k2, which according to the plot is approximately k = 75.

To begin with use another name for your procedure, e.g. MyPrime. The reason is that Prime (and prime) are protected:
type(Prime,protected);

Secondly, use set braces: { } NOT parentheses () for sets.

Thirdly, You need do give the output at the end.

Fourthly, I don't understand your line about r, so I cannot tell you what to do, but I can tell you that in any case you cannot refer to p[i] as that would mean that you will be referring to p[1], p[2],..., p[10] if the call MyPrime(1,10) is made.
So here the changes are made, but the line about r is removed (since I don't know your intention, so cannot repair):

restart;
MyPrime:= proc (a, b)
local i, p, r;
p := {};
#r:={};
for i from a to b do
  if isprime(i) then p := p union {i} end if;
end do;
p
end proc;

data:=[["a",13,"b",23],["c",2,"a",14],["d",4,"a",12],["a",8,"e",5],["f",3,"a",2]];
f:=L->if member("a",L,'p') then ["a",L[p+1],L[1..p-1][],L[p+2..-1][]] else L end if;
map(f,data);

The main problem is that M is a sequence with 5 elements, but your loops use M[n+1] = M[6].
Here I simply made M longer. You can adjust to obtain what it is you want.

f:=x->2^x;
n:=5;
M:=-2,-1,0,1,2,3;
P:=1;
for k to n+1 do
  L:=1;
  for i to n+1 do
    if i<>k then L:=L*(x-M[i])/(M[k]-M[i]) end if
  end do;
  P:=P+f(M[k])*L
end do;

I changed your code for s1 and s2 to use an inert integral (Int, not int). Also I used definite integrals. Surely the symbolic integration will fail (correction below) so numerical integration has to be used (just notice the RootOf coming out of solve).

s1 := Int(sqrt(1+(diff(subs(C, z1), x))^2), x=0..xP): #definite inert integral
s2 := Int(sqrt(1+(diff(subs(C, z2), x))^2), x=l..xP): # Ditto
s:=s1-s2:
eq1 := l0-s:
l0 := 40:
plot(eq1, H=10 .. 11, colour = red);
H := fsolve(eq1 = 0, H = 10.7 .. 10.8);
                               H := 10.7313620141804

##Plotting will be faster if you force plot to use fewer points:

plot(eq1, H=10 .. 11, colour = red,adaptive=false,numpoints=10);

######################Integration is not the problem:
restart;
z1:=H/mg*cosh(mg/H*(x+C11))-C21;
z2:=H/mg*cosh(mg/H*(x+C12))-C22;
tg_alpha1:=diff(z1,x);
tg_alpha2:=diff(z2,x);
r1:=H*subs(x=xP,tg_alpha2)-H*subs(x=xP,tg_alpha1)-P;
r2:=subs(x=0,z1);
r3:=subs(x=l,z2)-h;
r4:=subs(x=xP,z1)-subs(x=xP,z2);
mg:=2:l:=20:h:=5:xP:=l/5:P:=10:
simplify(Int(sqrt(1+tg_alpha1^2),x=0..xP)) assuming positive;
s1:=value(%);
simplify(Int(sqrt(1+tg_alpha2^2),x=l..xP)) assuming positive;
s2:=value(%);
C := solve({r1, r2, r3, r4}, {C11, C12, C21, C22}):
s:=eval(s1-s2,C):
eq1 := l0-s:
l0 := 40:
plot(eq1, H=10 .. 11, colour = red,adaptive=false,numpoints=10);

fsolve(eq1 = 0, H = 10.7 ); #Fast
#H := fsolve(eq1 = 0, H = 10.7 .. 10.8); #Takes a long time for some reason

It appears that there is a bug in int in Maple2015.1:

I5:=Int((x^2)/((x^4-2*x0*x^2+x0^2+1)^(3/2)),x=0..infinity);
I5a:=Student:-Precalculus:-CompleteSquare(I5,x0); #Just for convenience
value(I5a);
value(I5);

In both cases we get signum( some expression )*infinity.

The integral is clearly convergent for any real x0 since its integrand behaves as 1/x^4 for large x.

# You would probably be forced to use fsolve in any case, so you could do:

restart;
I5:=Int((x^2)/((x^4-2*x0*x^2+x0^2+1)^(3/2)),x=0..infinity);
I6:=(1/2)*Int(1/(x^4-2*x0*x^2+x0^2+1)^(1/2),x=0..infinity);
A0:=((-4/Pi)*((x0*I6-I5)/(x0*I5+I6)^(1/3)));
p:=proc(x00) evalf(eval(A0,x0=x00)) end proc;
p(1);
fsolve(A0=1,x0); #f=1 as an example
fsolve(A0=x0,x0); #f=x0 as another

PS. The bug is also present in Maple 15.01.
PPS. Maple12.02 behaves better.
value(I5) assuming x0>0;
outputs limit(expr,x=infinity), where expr is an expression in x and x0.

You want to find nontrivial solutions to a linear and homogeneous boundary value problem. Obviously there is the trivial solution, which exists no matter what 'a' is.
So you want to determine 'a' so that the problem has a nontrivial solution.
#The first version had the wrong eq4.
restart;
ode := diff(Y(x), x$4) = a^4*Y(x);
ics := Y(0) = 0, D(Y)(0) = 0, Y(l) = 0, (D@@2)(Y)(l) = 0;
res:=dsolve(ode);
eq1:=eval(rhs(res),x=0)=0;
eq2:=eval(rhs(res),x=l)=0
eq3:=eval(diff(rhs(res),x),x=0)=0;
eq4:=eval(diff(rhs(res),x,x),x=l)=0;#Edited
#You need to determine 'a' so that the linear system of homogeneous equations eq1,..,eq4 has a nontrivial solution:
A,z:=LinearAlgebra:-GenerateMatrix([eq1,eq2,eq3,eq4],[_C1,_C2,_C3,_C4]);
det:=simplify(LinearAlgebra:-Determinant(A)); #Needs to be zero
#Assuming that a = 0 is uninteresting:
det1:=op(3,det);
#det1 = 0 can be written as det2=0
det2:=cosh(a*l)*sin(a*l)-sinh(a*l)*cos(a*l);
plots:-logplot(eval(det2,l=2),a=0..10);
### Now we find the lowest 'eigenvalue' (i.e. a) when l = 2:
a2:=fsolve(eval(det2,l=2),a=3);
A2:=eval(A,{a=a2,l=2});
LinearAlgebra:-Rank(A2); #Returns 3 as it ought to.
LinearAlgebra:-GaussianElimination(A2);
fnormal~(%); #Removal of roundoff
LinearAlgebra:-LinearSolve(%,z,free=t);
eval(res,[_C1,_C2,_C3,_C4]=~convert(%,list));
Y2:=eval(rhs(%),{t[1]=1,a=a2});
plot(Y2,x=0..2);

###The next one:
a22:=fsolve(eval(det2,l=2),a=4);
A22:=eval(A,{a=a22,l=2});
LinearAlgebra:-Rank(A22); #Returns 3 as it ought to.
LinearAlgebra:-GaussianElimination(A22);
fnormal~(%);
LinearAlgebra:-LinearSolve(%,z,free=t);
eval(res,[_C1,_C2,_C3,_C4]=~convert(%,list));
Y22:=eval(rhs(%),{t[1]=1,a=a22});
plot(Y22,x=0..2);

Here is a way to find a solution for v(x) which continues after v(x) and v'(x) have become zero.

eq1:=op(solve(SYS[1],{u(x)}));
subs(eq1,SYS[2]);
eq2:=collect(%,diff);
res:=dsolve({eq2, v(-2) = 4, D(v)(-2) = 0,b(-2)=0. },numeric,discrete_variables=[b(x)::float],events=[[v(x)=0,b(x)=x]],output=listprocedure);
plots:-odeplot(res,[x,v(x)],-5..0);
V,V1,B:=op(subs(res,[v(x),diff(v(x),x),b(x)]));
x0:=B(-.1);
V(x0),V1(x0); #Nearly zero
#Now how does u behave if it is determined from eq1?
U:=proc(xx) eval[recurse](rhs(eq1),{v(x)=V(xx),diff(v(x),x)=V1(xx),x=xx}) end proc;
##U is discontinuous at x0, however:
U(x0);
U(x0-(1e-6));
plot(U,-5..0);

Your system SYS is an interesting example of a system with nonuniqueness. As pointed out in my comment above, when v(x) becomes 0 at x=x0 then v can be continued as zero and u can be continued as anything differentiable for x>x0.

factor(evalf(rho_poly)); #evalf
nops(%); #33
degree(rho_poly,rho); #32
So the answer is yes.

In both solve commands change the derivative of f4(x) to diff(f4(x), x$3) as that is the highest order derivative of f4 appearing in sys as well as in newsys.
I recognize your code as my code given to (I presume) a fellow student of yours for a similar system. You need to understand why you change from sys to newsys.

#### Note (edited)
It appears that you need to differentiate sys[4] also in order to solve for the highest derivatives:

newsys := {sys[1], sys[3], sys[5], sys[6], diff(sys[2], x), diff(sys[4], x)};

Obviously, the definition of bcs2 should contain the undifferentiated versions of the ones differentiated above.

I don't think that the two equations determine functions ya(x), yc(x).
I'm assuming that E and f are constants (and that by pi you mean Pi).

restart;
eq1:=(1-Pi*x/2)*(f/2/E)*yc-1/3=Int(t^4/sqrt((t^2-(f/2/E)*yc)^2+ya^2),t=0..1);
eq2:=x^2*yc*(1-Pi*x/2)-1/3=Int(t^4/sqrt((t^2-x^2*yc)^2+ya^2),t=0..1);
subs(yc=2*E/f*y1,eq1);
subs(yc=y2/x^2,eq2);
eq:=subs(y2=y,%);
##
Suppose for a given x there is a solution (ya,yc) to the system {eq1, eq2}.
Then y1 = yc*f/(2*E) and y2=yc*x^2 both satisfy eq above.
SUPPOSE that actually eq has a unique solution for y for given x and ya. If that is so then y1=y2, i.e.  yc*f/(2*E) =yc*x^2. Now yc=0 does not satisfy eq1 or eq2. Thus f/(2*E) =x^2, contradicting the assumption that f and E are constants (assuming of course that you want to do this for more than just one x).

I shall only give graphical evidence of unique solution to eq for given x and ya:
We may restrict ourselves to ya>=0.
Notice further that any solution y to eq is positive if x < 2/Pi and negative if x > 2/Pi.
#First a couple of test plots:
plots:-implicitplot(eval(eq,x=0),ya=0..4,y=0..2,gridrefine=2);
plots:-implicitplot(eval(eq,x=1),ya=-2..2,y=-2..2,gridrefine=2);
#Animations. You can adjust as you please, the number of frames in particular.
use plots in
  animate(implicitplot,[eq,ya=0..2,y=0..2],x=-3..2/Pi-.1,frames=15)
end use;
use plots in
  animate(implicitplot,[eq,ya=0..2,y=-4..0],x=2/Pi+.1..3,frames=5)
end use;

#It appears that for each given x and ya there is a uniquely given yc.

Is this what you are trying to do?

dsol2:=`union`(dsol[]);
eval(eq2,dsol2);

Not surprising that you get diff(s(t),t).
The same is obtained by
odetest(dsol,eq2);
which basically evaluates either eq2 (as here) or the difference between the left and right hand sides of the ode as in this case
odetest(dsol,eq2=diff(s(t),t));

Of course for eq3 and eq4 you get the same as for eq2:

eval([eq2,eq3,eq4],dsol2);
odetest(dsol,[eq2,eq3,eq4]);
odetest(dsol,[eq2,eq3,eq4]=~diff(s(t),t));

When you do

dsolve({Teq,UEQ});

you get a result that gives you 3 different solutions for T(eta):

1. T(eta) = 0, which you cannot use because of your boundary conditions
2. T(eta) = _C3*eta+_C4
3. A complicated looking expression giving T(eta) in terms of a RootOf. Out of this one you are surely not going to get an analytic answer.

u(eta) is given as an unevaluated double integral involving T(eta).

So the best (or only) hope is for T(eta) = _C3*eta+_C4 to give us an expression for u(eta) where the integrals can be computed.
By using the boundary conditions we get T(eta) = eta.
After that we try

dsolve(eval(UEQ,T(eta)=eta));

We notice that we no longer have a double integral, but an integral, which Maple cannot do as it stands:
value(%);

You can try
dsolve({eval(UEQ,T(eta)=eta),D(u)(1)=0,u(0)=L*D(u)(0)}); #You don't get anything
dsolve({eval(UEQ,T(eta)=eta),u(0)=L*D(u)(0)}); #You do get something but an integral is still there
dsolve({eval(UEQ,T(eta)=eta),D(u)(1)=0}); #Also here you get an integral
############
## Numerical solution:
gama1:=0.2:
rhop:=5180:
rhobf:=998.2:
a[mu1]:=39.11:
b[mu1]:=533.9:
a[k1]:=7.47:
b[k1]:=0:
L:=1: N_bt:=1: #Not supplied by you.
res:=dsolve({Teq,UEQ,D(u)(1)=0,u(0)=L*D(u)(0),T(1)=1,T(0)=0},numeric);
plots:-odeplot(res,[eta,T(eta)-eta],0..1);
## In this case we do get a solution where T(eta) = eta. See note below ******
#This gives us not much of a hint as to the form of the value of the integral for u:
plots:-odeplot(res,[eta,u(eta)],0..1);
res(0);
res(1);
####
***Note: This is not surprising in view of the result of
`dsolve/numeric/bvp/convertsys`({Teq,UEQ},{D(u)(1)=0,u(0)=L*D(u)(0),T(1)=1,T(0)=0},[T,u],eta,l);
Or the result of:
solve(Teq,diff(T(eta),eta,eta));