@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);

@minhthien2016  
simplify((1/2)*ln(1+tan(x)^2))  assuming cos(x)>0;
                                -ln(cos(x))

@minhthien2016 

restart; 
int(sin(x)^n, x) assuming n::posint:
eval(%, n = 7); 
expand(simplify(value(%)));
algsubs(cos(x)^2=1-sin(x)^2, %);
factor(%);

@student_md 
Tm2 is present in your new code, but it is not defined. From your previous version it seems that  Tm2  is the same as Tm. Correct this.

The first error in your code occurs in the following line:

coz:=fsolve(sys);
indets(sys, name);

Error, (in fsolve) number of equations, 16, does not match number of variables, 21
{P_nn3, c[1, 1], c[1, 2], c[1, 3], c[1, 4], c[2, 1], c[2, 2],  c[2, 3], c[2, 4], c[3, 1], c[3, 2], c[3, 3], c[3, 4], [4, 1],  c[4, 2], c[4, 3], c[4, 4], t[1], t[2], t[3], t[4]}
 

@Jalale  See the corrected file:

QuestionSimulation_new1.mw 

@gaurav_rs  Maple often does not take into account restrictions on variables imposed by  assume  command.

An example:

assume(x>Pi/2, x<Pi):
solve(sin(x)=1/2);
            Warning, solve may be ignoring assumptions on the input variables.
                              1   
                              - Pi
                              6   
 

@Magma  See help on  ?remember . There everything is explained in detail with examples.

@vv   ``( ... )  construction prevents the automatic opening of the parentheses.

@mehdi jafari   Probably no.

@digerdiga 

f:=ln((1-x)^2*(x+1)^2/((-I*x-I+sqrt(-x^2+1))^2*(I*x+I+sqrt(-x^2+1))^2))*ln((1-x)^2*(x+1)^2/((-I*x-I+sqrt(-x^2+1))^2*(I*x+I+sqrt(-x^2+1))^2))*ln((1-x)^2*(x+1)^2/((-I*x-I+sqrt(-x^2+1))^2*(I*x+I+sqrt(-x^2+1))^2))*ln((1-x)^2*(x+1)^2/((-I*x-I+sqrt(-x^2+1))^2*(I*x+I+sqrt(-x^2+1))^2)) + ln((1-x)^2*(x+1)^2/((-I*x-I+sqrt(-x^2+1))^2*(I*x+I+sqrt(-x^2+1))^2))*ln((1-x)^2*(x+1)^2/((-I*x-I+sqrt(-x^2+1))^2*(I*x+I+sqrt(-x^2+1))^2)) + ln((1-x)^2*(x+1)^2/((-I*x-I+sqrt(-x^2+1))^2*(I*x+I+sqrt(-x^2+1))^2)) + ln((1-x)^2*(x+1)^2/((-I*x-I+sqrt(-x^2+1))^2*(I*x+I+sqrt(-x^2+1))^2))*ln((1-x)^2*(x+1)^2/((-I*x-I+sqrt(-x^2+1))^2*(I*x+I+sqrt(-x^2+1))^2));
subsindets(f, `function`, p->applyop(t->numer(t)/expand(denom(t)),1, p));
simplify(%);
expand(%);