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

@ecterrab  But this is not an automatic solution. In fact, you manually wrote an expansion, which I hoped to get from Maple.

@nm  Unfortunately, I was not able to verify by differentiating the result in my answer above for an arbitrary  n  and get  sin(x)^n . However, the correctness of this result is beyond doubt. Here is a simple check for 100 values of  n :

restart;
R:=int(sin(x)^n, x) assuming n::posint;
seq(is(value(eval(R, n=k))=int(sin(x)^k, x)), k=1..100);