@vv  In fact, as a formal symbolic result, we get the same thing. But if we try to assign some numeric constants to the symbols, then of course we get an error:

Ec := (Ems+I*Eml)*(1+((Ems+I*Eml)/Ef-1)*Zeta*phi/((Ems+I*Eml)/Ef+Zeta))/(1-((Ems+I*Eml)/Ef-1)*phi/((Ems+I*Eml)/Ef+Zeta)):
a:=simplify(Re(Ec)) assuming positive;
b:=simplify(Im(Ec)) assuming positive;
EC:=eval(Ec,Zeta=ZETA):
simplify(evalc([Re,Im](EC)));
eval(%, ZETA=Zeta);
is(a=%[1] and b=%[2]);
Zeta:=1;

@bsoudmand  This is the output of my code (in Maple 2018.2):

In the text form:

a:=(-Ems*Zeta*(Ef^2-2*Ef*Ems+Eml^2+Ems^2)*phi^2+((Zeta-1)*Ems^3+Ef*(Zeta-1)^2*Ems^2-(Zeta-1)*(Ef^2*Zeta-Eml^2)*Ems-Ef*Eml^2*(Zeta+1)^2)*phi+Ems*(Ef^2*Zeta^2+2*Ef*Ems*Zeta+Eml^2+Ems^2))/((Ef^2-2*Ef*Ems+Eml^2+Ems^2)*phi^2+(-2*Ems^2-2*Ef*(Zeta-1)*Ems+2*Ef^2*Zeta-2*Eml^2)*phi+Ef^2*Zeta^2+2*Ef*Ems*Zeta+Eml^2+Ems^2);

b:=-Eml*(Zeta*(Ef^2-2*Ef*Ems+Eml^2+Ems^2)*phi^2+(Ef*(Ef-2*Ems)*Zeta^2+(-Ef^2-Eml^2-Ems^2)*Zeta-2*Ef*Ems+Eml^2+Ems^2)*phi-Ef^2*Zeta^2-2*Ef*Ems*Zeta-Eml^2-Ems^2)/((Ef^2-2*Ef*Ems+Eml^2+Ems^2)*phi^2+(2*Ef*(-Ems+Ef)*Zeta+2*Ef*Ems-2*Eml^2-2*Ems^2)*phi+Ef^2*Zeta^2+2*Ef*Ems*Zeta+Eml^2+Ems^2);

@vv 

SymFun(ln((x-1)/(x+1)), x);
plot(ln((x-1)/(x+1)), x=-4..4, -4..4);

SymFun(sqrt(x^2), x);
plot(sqrt(x^2), x=-2..2);

SymFun(abs(x)^(1/2), x);
plot(abs(x)^(1/2), x=-1..1);
 

@emendes  You can use  makeproc  option. See help on  rsolve  command for this.

@Lali_miani  If you enter in 2d math mode any of  i, j, I  from the Comman Symbols palette, it will work as the imaginary unit. But from the keyboard you have to enter the imaginary unit as  I . I almost never use palettes and only work from the keyboard in 1d math mode (Maple input). 

@minhthien2016  You can make arrows at the ends of the axes of coordinates and labels in the same place if you use  plots:-arrow  and  plots:-textplot  commands.

Example:

Arrow_x:=plots:-arrow([2.3,0],[0.15,0], width=0, head_width=0.12, head_length=0.15):
Arrow_y:=plots:-arrow([0,4.3],[0,0.15], width=0, head_width=0.12, head_length=0.15):
Labels:=plots:-textplot([[2.4,-0.2,"x"],[-0.2,4.4,"y"]], font=[times,16]):
P:=plot(x^2, x=-2..2, color=red, thickness=2, labels=["",""]):
plots:-display(Arrow_x,Arrow_y,P,Labels);

Output:
                         

Edit.

@bliengme 

Add the line

convert(Q, exp);

@EB1000 The  DirectSearch  package is not included in Maple and must be downloaded from the Maple Application Center from  here

@tomleslie  Thanks for this. In my opinion it is the best solution.

@radaar  I did not understand the meaning of what you posted. But in general, I think that using loops is more efficient, but you need to increas  Digits  (for example take Digits:=15)  to compensate for the loss of accuracy.

@radaar Note that double  add  requires more time and memory, but the calculation itself is more accurate. In order to achieve the same accuracy with double  for-loop  you have to increase  Digits , but then both time and memory will increase. You need to check it all out. 

@mmcdara  I do not think your solution is shorter. You have written a special procedure for this, which uses a number of commands from  Statistics  package. I simply use  Composition  procedure, written several years ago and not directly related to this problem. The solution itself takes 1 line of the code and is about 50 times faster (see below):
 

Download SumOfDice_new.mw

@torabi  You have already been offered to use  DirectSearch  package. It is not part of Maple and must be downloaded from Maple Application Center from  here

@mehran rajabi 
Replace the line

Sys:={seq(u[k]=1+add((x[k]*x[i-1]*u[i-1]+x[k]*x[i]*u[i])/2*h, i=1..N), k=1..N)}:

with the lines

f:=i->x[k]*x[i]*u[i]:
Sys1:={seq(u[k]=1+h/3*(f(0)+4*add(f(2*i-1), i=1..N/2)+2*add(f(2*i-2),i=2..N/2)+f(N)), k=1..N)}:

Thank you all for the answers. Apparently there is no simple solution and it is easier to simply record this expansion handly.