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;

@DavidABurton Yes, I can confirm that on Maple 2015.2 and using Windows 10 the plot looks as you show.
On the same machine I tried Maple 18 and Maple 2016.2. In both cases the plots looked much better than the Maple 2015 version, but the Maple 2016 plot didn't look quite as nice (in some sense) as the Maple 18 version.
I also looked at an earlier version, Maple 15. No problem there either.

Maybe one of the experts in these matters can tell us what is going on in Maple 2015.2 w.r.t. graphics.
I suspect that it is a "well known" problem.
 

Since you have tried several things you should upload the worksheet so somebody can look at it.
Use the fat green arrow in the editor.

@nm Here is an non-numeric version:

restart;
ode:=diff(y(x),x$2)+lambda(x)*y(x)=0,diff(lambda(x),x)=0;
bc:=y(0)=0,y(L)=0,D(y)(0)=y1;
res:=dsolve({ode,bc});

Since _2*Z1+_B1 can take on all integer values you can continue with:

subs(2*_Z1+_B1=n,res);
simplify(%) assuming L>0,n>0;

where I have also used yhat we may as well assume that n > 0.