You need x->sqrt(-x) and x->-sqrt(-x).

I only looked at the first integral, which I denoted A, but in the second integral x and z are just replaced by
x+l*cos(alpha) and z+l*sin(alpha), respectively.
Maple's answer is using piecewise (see the help for piecewise).
Unfortunately, the answer is not quite right as we shall see.
restart;
A:=Int(1/sqrt((x-R*cos(theta))^2+(y-R*sin(theta))^2+z^2), theta = 0 .. 2*Pi);
#First I try with actual values for R, x, y, and z. This will be a test case.
value(eval(A,{R=1,x=2,y=3,z=4})); #The exact result
evalf(%); #The decimal version of the exact result
evalf(eval(A,{R=1,x=2,y=3,z=4})); #Numerical integration instead
#As the last two results agree we may believe the exact result in this test case.
#Now trying in general:
value(A); ## Result = 0. Definitely wrong!
#Now we try first simplifying as you did:
simplify(A);
res:=value(%): #This the the piecewise result, which says that under some conditions the result is undefined and otherwise it is given by a certain expression (which I shall call U) involving elliptic functions.
#We test the result on our test case:
eval(res,{R=1,x=2,y=3,z=4}); # Result undefined, not good!
#The 'otherwise' part of the piecewise result we can get as the last operand in res:
U:=op(-1,res):
eval(U,{R=1,x=2,y=3,z=4});
evalf(%);
#Except for roundoff error this actually is correct.
So maybe U is the correct answer at least for real values of x,y,z,R (not all zero).
A neater version of U can be found this way:

simplify(U) assuming real,R>0;
algsubs(x^2+y^2=r^2,%); #Using r = sqrt(x^2+y^2)
U1:=simplify(%) assuming r>0;
#And more work will make it even simpler:

u1:=R^2+2*R*r+r^2+z^2;
u2:=R^2-2*R*r+r^2+z^2;
subs(u1=p,u2=m,U1);
simplify(%) assuming positive;
U2:=subs(p=(R+r)^2+z^2,m=(R-r)^2+z^2,%);

I have included the empty set in the output:

p:=proc(S::set) local T1,T2;
  if nops(S)=1 then
    S union {{}}
  else
    {seq(seq(seq(`union`(T1,T2),T1=procname(S[2..-1])),T2=S[i..i]),i=1..nops(S)),{}}
  end if
end proc;
p({A,B,C,E});
p({{a},{b},{c,b,67},{8}});

Changing the time scale may be what you need:

inits:={a[1,1](0)=0.1,D(a[1,1])(0)=0.1,a[2,1](0)=0.1,D(a[2,1])(0)=0.1,a[1,2](0)=0.1,D(a[1,2])(0)=0.1,a[2,2](0)=0.1,D(a[2,2])(0)=0.1,a[1,3](0)=0.1,D(a[1,3])(0)=0.1,a[2,3](0)=0.1,D(a[2,3])(0)=0.1,a[3,3](0)=0.1,D(a[3,3])(0)=0.1,a[3,2](0)=0.1,D(a[3,2])(0)=0.1,a[3,1](0)=0.1,D(a[3,1])(0)=0.1};
sys:={seq(cat(eq,i),i=1..9)}:
SYS:=solve(sys,{seq(seq(diff(a[i,j](t),t,t),i=1..3),j=1..3)}):
seq(eval(subs(t=0,convert(rhs(SYS[i]),D)),inits),i=1..9);
#Choosing the new time s=10^6*t:
SSYS:=PDEtools:-dchange({t=10^(-6)*s},SYS):
S,R:=selectremove(has,inits,D);
map(c->lhs(c)=rhs(c)*10^(-6),S);
Sinits:=% union R;
MMM := dsolve(SSYS union Sinits,numeric);
MMM(10.234);
plots:-odeplot(MMM,[s,a[1,1](s)],0..25);
plots:-odeplot(MMM,[s,a[1,2](s)],0..25);
plots:-odeplot(MMM,[s,a[1,3](s)],0..25);
plots:-odeplot(MMM,[s,a[2,1](s)],0..25);
PLTS:=seq(seq(plots:-odeplot(MMM,[s,a[i,j](s)],0..25,caption=a[i,j]),j=1..3),i=1..3):
plots:-display(PLTS,insequence=true,labels=[s,a]):

If you look at the output from the assignment to d you will notice that the version of x[E] and y[E] you have there are not evaluated although you have assigned to `x__E` and `y__E`. These are not the names used in the d-assignment.
Once that issue is resolved you may still not have luck with fsolve, but it will not give you an error.

Of course it is easily checked that the expressions are equal. It is much harder to make Maple go from the first to the second.
Student:-Precalculus:-CompleteSquare is made for that kind of thing, but doesn't help all that much.

restart;
u:=(2*delta^2-2*beta*delta+alpha^2-2*alpha*delta+beta^2)*N^2+2*gamma*(beta+alpha-2*delta)*N+2*gamma^2;
Student:-Precalculus:-CompleteSquare(coeff(u,N^2),[alpha,beta]);
#So we do this first:
collect(u,N,x->Student:-Precalculus:-CompleteSquare(x,[alpha,beta]));
#Introduce temporary variables:
subs(alpha=delta+ad,beta=delta+bd,%);
Student:-Precalculus:-CompleteSquare(%,[ad,bd]);
simplify(%,power);
subs(ad=alpha-delta,bd=beta-delta,%);

Here are 3 ways:

res:=solve({x-y = -1, 2*x-y = 3}, {x, y});
#First way
u:=subs(res,[x,y]);
u[1];
u[2];
#Second way and using capital X, Y instead of x,y:
X,Y:=op(subs(res,[x,y]));
X;
Y;
#The third way assign to x and y, which I personally almost never would do:
assign(res);
x,y;

Doesn't this work for you:

Z:=[seq(evalf(sin(k)^2),k=1..12)]; #An example of a list of reals
sum(ln(1-exp(Z[i]/T)),i=1..12); #Using sum
#But for concrete summations i.e. with known bounds, here i = 1..12 'add' is the recommended one:
add(ln(1-exp(Z[i]/T)),i=1..12);

In your first example you can include maxfun=0 (or set a high number) as in
sol:=dsolve({eqbz,ic},numeric,output=listprocedure,parameters=[l],maxfun=0); #No range

I don't know what you mean with these lines:

ysol(parameters=[l=x]):

current:=convert(ysol(parameters,`global`):

For the other example you can use the procedural version, see ?dsolve,numeric,ivp
You cannot use the parameters option in that case, but you can introduce the parameter l as a second function satisfying diff(l(x),x)=0, l(x0)=the l-value desired.
Something like:

solproc:=proc(N,x,Y,YP) local yeq,p;
    yeq:=x^3*evalf(Int(p^2/(exp(x*sqrt(p^2+1))+1),p=0.001..10^4));
    YP[1]:=-Y[2]/x^2*(Y[1]^2-yeq^2);
    YP[2]:=0
end proc:
#With l = 2.3:
sol:=dsolve(numeric,number=2,procedure=solproc,initial=Array([1,2.3]),start=0.01,
   output=listprocedure,procvars=[y(x),l(x)]);

Alternatively, you could do:

restart;
solproc:=proc(N,x,Y,YP) local yeq,q;
    yeq:=x^3*evalf(Int(p^2/(exp(x*sqrt(p^2+1))+1),p=0.001..10^4));
    YP[1]:=-L/x^2*(Y[1]^2-yeq^2);
end proc:
sol:=l->dsolve(numeric,number=1,procedure=subs(L=l,eval(solproc)),initial=Array([1]),start=0.01,
   output=listprocedure,procvars=[y(x)]);
res:=sol(2.3); #With l = 2.3


There are numerical problems however, so there is more work to be done or options to insert into dsolve.

Leave out the 'do' following right after 'else' and consequently leave out one 'od'.

if 3>4 then
   print(blah)
else
   i:=0;
   while i < 3 do
      i:=i+1;
      print(i)
   end do
end if:

Or use a for loop:
if 3>4 then
   print(blah)
else
   for i from 0 to 20 while i^2 < 8 do
      print(i);
      print(i^2)
   end do
end if:

By dividing the matrix B with the matrix A do you mean solving the equation A.X = B or du you mean solving X.A = B?
In the first case X = A^(-1).B in the second X = B.A^(-1) both assuming that A is invertible (regular).
This would be Maple syntax too.

I guess you don't mean elementwise division, but just in case you do this can be done like in this example:

A:=Matrix(2,symbol=a);
B:=Matrix(2,symbol=b);
A/~B;

The only meaningful equilibrium point is not (necessarily) number 3 as you have written:
fix:=subs({theta(t)=theta, v(t)=v},%[3]);
In fact on my machine (and with Maple 17) it came out as number 1, so the line should be:
fix:=subs({theta(t)=theta, v(t)=v},%[1]);
#With that kind of code you have to be very careful: Does the order remain consistently the same or not?
#You could do this instead:
solve( [rhs(phug[1]), rhs(phug[2])=0,v(t)>0], [theta(t),v(t)]);
#Now you have 2 "solutions", but the speeds are the same and the angles differ by Pi, but the one is in fact wrong as a simple check will reveal.
The correct one is
[theta = -.3805063771, v = .9635749534]

You could use Maximize from the Optimization package instead of maximize. F is just a simple factor in your expression. Just leave it out to begin with.

evalf(-N/(b*h)+((-EI*(diff(u(x), x, x)))*.5)*h/(0.10e-10*Iice));
collect(%,F);
plot(op(1,%) ,x=0..l);
Optimization:-Maximize(eval(-N/(b*h)+((-EI*(diff(u(x), x, x)))*.5)*h/(0.10e-10*Iice), F = 1), x = 0 .. l);

p:=plot(---what you got ---);

tf:=plottools:-transform((x,y)->[y,x]):
tf(p);
plots:-display(tf(p),labels=["u(y)","y"]);

You have got a wrong sign in your last equation.
From your equations it appears that the 3 masses are 10, 5, 10 and the 4 spring constants each have the value 5. Thus the system will be:

a1 := {10*(diff(diff(a(t), t), t)) = -5*a(t)+5*(b(t)-a(t)), 5*(diff(diff(b(t), t), t)) = -5*(b(t)-a(t))+5*(c(t)-b(t)), 10*(diff(diff(c(t), t), t)) = -5*(c(t)-b(t))-5*c(t)};
This can be written in the form M.X'' + K.X = 0 or equivalently X'' + M^(-1).K.X = 0, where X = <a,b,c> and where M is the diagonal matrix containing the 3 masses and K the matrix
Matrix([[10, -5, 0], [-5, 10, -5], [0, -5, 10]])
The eigenvalues of M^(-1).K are 1 , (3+sqrt(5))/2, (3-sqrt(5))/2 thus positive as they ought to be.
The solution to your problem can be done directly by dsolve like this:

res:=dsolve(a1 union {a(0) = 0, D(a)(0) = 1, b(0) = 1, D(b)(0) = 2, c(0) = 3, D(c)(0) = 3});
and plotted like this:
plot(subs(res,[a(t),b(t),c(t)]),t=0..50);