I know of no word command for squaring. But you can easily make one yourself, if you feel the need.

L:=[7,9,13,evalf(Pi),a,b];
sq:=x->x^2:
sq~(L);
map(sq,L); #Alternative
##And without a name:
map(x->x^2,L);
(x->x^2)~(L); #alternative

## Notice that sqrt is not just defined as sqrt:= x-> x^(1/2); as will be seen here:
showstat(sqrt);

## And a simple example of the difference:
4^(1/2);
sqrt(4);

Changing the order of integration in the double integral below is done by hand.
(Revised version)
restart;
intE := diff(u(x), x) = 1+x-x^2+int((x-t)*u(t), t = 0 .. x);
init:= u(0) = 3;
map(int,intE,x=0..y,continuous); #'continuous' necessary for left hand side only
eval(%,{init,continuous=NULL});
combine(IntegrationTools:-Expand(%));
subsop([2,1]=int(int(u(t)*x-u(t)*t, x = t .. y), t = 0 .. y), %); #Changing order of integration
factor(%);
eq:=isolate(subs(y=x,%),u(x));
diff(eq,x); #Check
eval(eq,x=0); #Check

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.