Could you provide us with all the code in the form of text or an uploaded worksheet?
Images don't help much.
Nobody is going to type the stuff and they won't even know how temp2 was defined or what happened before.

@ernilesh80 Since you don't know how large the two arrays have to be you can just define them as arrays:
pd_arr:=Array();
id_arr:=Array();
###
Then in the assignments to pd_arr and id_arr in the loops you must use parentheses instead of square brackets as in
pd_arr(pd, 1) := op_value[1]; pd_arr(pd, 2) := op_value[2]; pd_arr(pd, 3) := op_value[3]; pd_arr(pd, 4) := op_value[4];
Do similarly for id_arr.
The loops take a while since
nops~([T_arr,E_arr,W_arr,p_arr]);
`*`(op(%)); 
returns 7436529.
I stopped the computation after a while and inspected both arrays by just doing:
pd_array;  # at that time a 1..8667x1..4 array
id_array ;  # at that time a 1..46980x1..4 array, i.e. larger than the 1..10000x1..4 you used.
## Incidentally, it may not be more efficient to do this, but it is shorter:
pd_arr(pd, 1..4) := Array(op_value);
and similarly for id_arr.

@memdream You should be aware that (-1)^(1/4) in Maple is the principal root:
evalc( (-1)^(1/4) );
                                 (1/2)*sqrt(2)+(1/2*I)*sqrt(2)

@Rouben Rostamian  Yes, when I paste the 2D input on a 1D input line I get
exp(4.605170186*`9`)
The factor `9` is understood by Maple as a name, not as the integer 9.

Since nobody is going to type your expression from the image you should upload a worksheet using the fat green arrow in the MaplePrimes editor.

Why not give us the specific function? (As text or an uploaded worksheet!)

@Rouben Rostamian  Yes, you are right. That is why the -1 is repeated, since res still has the value from i = -1.

@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?