@tomleslie gamma is Euler's constant, thus gamma is approximately 0.5772156649.

If x and y are not assigned values of any kind (using the assignment operator :=  as in x:=23 ) then the result of the statement
x+y;
is just that x+y. What else would you expect?
If you do
x:=a;
y:=78;
then
x+y;
evaluates to a+78 assuming that the variable a is unassigned.

For us to understand what you actually did you should upload a worksheet using the fat green arrow in the MaplePrimes editor.

@MrMarc To confirm that you need unevaluation quotes, try this:

rand()$10;
## and then
'rand()'$10;
## Take a look at the results from your Maple 2017.

The sequence operator `$` evaluates its arguments fully as the first thing (the normal behavior in Maple).
Thus rand() is evaluated once in the first case above and the result is a sequence of 10 identical numbers.
In the second case full evaluation also takes place, but full evaluation of an expression with unevaluation quotes just means removal of the quotes, thus when `$`gets to do its actual work it sees rand() and not a number as in the first case.

If you use seq instead you won't have this 'problem'. seq has special evaluation rules (as do e.g. add and mul).
With or without unevaluation quotes this works as intended:

seq(rand(),i=1..10);

I only have access to Maple 12 at this moment, but your test

time(rand()$1000000);

needs to be
time('rand()'$1000000);   # Unevaluation quotes

in Maple 12. Otherwise rand() executes only once and is then just repeated 1000000 times.
On my little Asus Eee PC with (I think) 2 GB RAM I get 37s for the latter command, but only 0.1s for the first!
 

@kambiz1199 The matrix equation you now have is actually worse than the previous equation. The 0 on the right hand side of course represents the 4x2 zero matrix. Without even requiring the unknowns to be real you now have 8 equations linear and nonhomogeneous equations with 7 unknowns. This system doesn't have any solution.

@MrYouMath Matrix form may mean different things to different people and to the same person different things at different occasions, but here is one version:
 

restart;
odeSys := {diff(x(t),t$2)+diff(x(t),t)+x(t)=f(t),diff(y(t),t$2)+2*diff(y(t),t)+3*y(t)=g(t)};
res:=DEtools[convertsys](odeSys,{x(t0)=x0,D(x)(t0)=x1,y(t0)=y0,D(y)(t0)=y1},[x(t),y(t)],t,Y,YP);
A,b:=LinearAlgebra:-GenerateMatrix( rhs~(res[1]), [seq(Y[i],i=1..4)]);
SYS:=diff(Y(t),t)=A.<seq(Y[i](t),i=1..4)> - b; 
ICS:= Y(res[3]) = <res[4]>;

Notice that if you leave out the optional names (here I used Y and YP) then local names (in fact also called Y and YP) will be used. That will mean that the subs command won't work as expected.

Just a comment about your loop:
You have set pr:=6.72;
Your loop reads
for j to nops(pr) do  (...content...)  end do;
That doesn't make much sense to me. pr is just the floating point number 6.72, so nops(pr) will be two because the operands of 6.72 are the mantissa (672) and the exponent (-2) of 6.72.
Try op(pr);
So forget about the loop and just try dsolve. I haven't had any luck.
I suppose you have a reason to believe that your boundary value problem has a solution?
If so you might also have an idea about the looks of such a solution which could be used as an approximate solution.

On your original expression normal and radnormal work fine:
 

u:=16*a^8*B/((dz*L*sqrt(s)*sqrt(s+c)*sqrt(L^2*s*(s+c)*dz^2+4*a^2)+L^2*s*(s+c)*dz^2+2*a^2)^2*(-dz*L*sqrt(s)*sqrt(s+c)*sqrt(L^2*s*(s+c)*dz^2+4*a^2)+L^2*s*(s+c)*dz^2+2*a^2)^2);
normal(u,expanded);
radnormal(u);
radnormal(u,rationalized);

On your new example I don't get anything that is shorter even with i = 1. There is an A in your larger example that is not present in your original. 
Why do you think that the new expression can be expressed without square roots?

Since we don't know what MODEL, VARS, DOMAIN, RANGE, [IC1, IC2] are you should upload a worksheet using the fat green arrow in the editor; alternatively give us as text the contents of  MODEL, VARS, DOMAIN, RANGE, [IC1, IC2] .

Executing your code after correcting a misspelling of MultiSeries (capital S) I got:

asympt(ii_inf,x,3);

Error, (in asympt) unable to compute series
MultiSeries:-asympt(ii_inf,x,1);
Error, (in MultiSeries:-multiseries) unable to sort exponents, {s, 2-s}

Nobody (I'm sure) would want to write the code seen in the image in your link.
So upload a worksheet using the fat green arrow in the MaplePrimes editor.

@Markiyan Hirnyk For your new example I get the exact same output for the two versions:

restart; 
int(ln(-a^2*x^2+1)/sqrt(-x^2+1), x = 0 .. 1);
restart;  
A := Int(ln(-a^2*x^2+1)/sqrt(-x^2+1), x = 0 .. 1):
value(A);

I think this is a limitation and not a bug.
The limitation is also shown in a separated version, where w = u-v and z = u+v. The two resulting pdes for w and z can be solved independently. The pde for w gives no problem, but the one for z shows the same problem as the given system.
 

restart;
sys:={diff(u(x, t), t)-diff(v(x, t), x)+u(x, t)+v(x, t) = (1+t)*x+(x-1)*t^2, diff(v(x, t), t)-diff(u(x, t), x)+u(x, t)+v(x, t) = (1+t)*x*t+(2*x-1)*t};
##
icbs:= {u(0, t) = 0, u(x, 0) = 0, v(0, t) = 0, v(x, 0) = 0};
## 
pde1:=eval(sys[1]-sys[2],u(x,t)=v(x,t)+w(x,t));
solw:=pdsolve(pde1,{w(x,0)=0,w(0,t)=0},numeric,time=t,timestep=0.01,spacestep=0.01,range=0..1);
solw:-plot3d(x = 0 .. 1, t = 0 .. 1); #No problem
pde2:=eval(sys[1]+sys[2],u(x,t)=-v(x,t)+z(x,t));
solz:=pdsolve(pde2,{z(x,0)=0,z(0,t)=0},numeric,time=t,range=0..1,timestep=0.01,spacestep=0.01);
solz:-plot3d( x = 0 ..1, t = 0 .. 1); #Problem

@Markiyan Hirnyk For an equal handling of int and value(Int(...)) clearly assuming should be applied when the integration kicks in.
Here is a more direct version:

restart;
int(ln(-a^2*x^2+1)/sqrt(-x^2+1), x = 0 .. 1) assuming a>0, a<1;
restart;
value(Int(ln(-a^2*x^2+1)/sqrt(-x^2+1), x = 0 .. 1)) assuming a>0, a<1;

@Markiyan Hirnyk Clearly Kitonum is quite right.
Your example certainly doesn't prove he is wrong.
The last part should have been

A := Int(ln(-a^2*x^2+1)/sqrt(-x^2+1), x = 0 .. 1);
res := value(A) assuming a>0, a<1;
which returns unevaluated as did this first one:
int(ln(-a^2*x^2+1)/sqrt(-x^2+1), x = 0 .. 1) assuming a>0, a<1;