with(plottools):
A:=rectangle([0, 0], [10, 10], style=line, thickness=2,color = red):
f:=p->rotate(p, arctan(1/3), [5,5]):
g:=p->homothety(p, sqrt(10)/4, [5,5]):
h:=g@f:
plots:-display(seq((h@@n)(A), n=0..10));

It looks like a bug.

3 different results. The plot confirms the correctness of the latter:

z:=exp(I*t):
f:=z^(1/2)*diff(z, t):
int(f, t=0..2*Pi);
int(evalc(f), t=0..2*Pi);
int(simplify(evalc(f)), t=0..2*Pi);
plot([Re(f), Im(f)], t=0..2*Pi, color=[red,blue], scaling=constrained);

Addition. For the second example, all methods (in Maple 2016.2) yield the same result:

f:=1/2+exp(I*t):
evalf(Int(f, t=0..2*Pi));
int(f, t=0..2*Pi);
int(evalc(f), t=0..2*Pi);
int(simplify(evalc(f)), t=0..2*Pi);
plot([Re(f), Im(f)], t=0..2*Pi, color=[red,blue], scaling=constrained);
 

Use  Array  and  Matrix  instead  array  and  matrix:

U := Array(1 .. 5, 1 .. 5);
convert(U, Matrix);
for i to 5 do
for j to 5 do
if i = j then U[i, j] := 1 else
U[i, j] := 0; 
end if end do end do: 
U;

A more compact ways for solution the same task:

U := Array(1 .. 5, 1 .. 5, (i, j) ->`if`(i = j, 1, 0));

or

U := Matrix(5,  (i, j)-> `if`(i = j, 1, 0));

Just assigning  Fh1  and  theta  values from the first solution:

sol := {Fh1 = 121.6477702, theta = 0.9606764638}, {Fh1 = -121.6477702, theta = -2.180916190}:
assign(sol[1]);
Fh1, theta; 
                                         121.6477702, 0.9606764638
 

If you want to save the obtained points  (x,y)  in a list to use later, for example, for plotting the function  y(x), you can assign the resulting values to indexed variables  x[i], y[i]  at each step of the loop:

restart;
x[0]:=0: y[0]:=1: xf:=1: n:=10:
h:=evalf((xf-x[0])/n):
f:=(x,y) -> x+y:
for i from 1 to n do 
    k:=f(x[i-1],y[i-1]);
    y[i]:=y[i-1]+h*k;
    x[i]:=x[i-1]+h;
od:
L:=[seq([x[i],y[i]], i=0..n)];   # List of received points
dsolve({diff(g(t),t)=t+g(t), g(0)=1}, g(t));  # The exact solution
assign(%):
plot([L, g(t)], t=0..1, color=[red, blue],  labels=["x", "y(x)"]);  # The plots of the approximate (red) and the exact (blue) solutions

Should be:

1.  exp(x)  instead of  e^x

2. A semicolon instead of a colon for displaying

3. A multiplication sign was missed in one place:

with(VectorCalculus):
with(LinearAlgebra):

Determinant(Wronskian([exp(2*x), exp(-x), x*exp(-x), x^2*exp(-x)], x));
                                          
                                         -54*exp(2*x)*(exp(-x))^3
 

NI.mw

First I simplified your function, because for  x>= 163  it is 0. This function is piecewise linear. Further, the integral is calculated as the sum of the areas of ordinary trapezoids. The calculation is verified by using  int  command:

f:=simplify(f);
L:=convert(f, list);
n:=nops(L);
L1:=[seq(L[2*k-1], k=1..(n-1)/2)];  # List of nodal points
P:=rhs~(L1);
F:=[seq(unapply(L[2*k],x), k=2..(n-1)/2)];
add((F[i](P[i])+F[i](P[i+1]))/2*(P[i+1]-P[i]), i=1..nops(F));  # Calculation of the integral by direct summation
int(f, x=0.594..163);  # Check by int command
plot(f , x=0.594..163);

Edit.

restart;
f:=(x,y)->x^2/4+y^2/8+1:
P:=[1,1,11/8]:
Surface:=plot3d(f(x,y), x=-1..3, y=-1..3, color="Yellow"):
Plane:=plot3d(D[1](f)(P[])*(x-P[1])+D[2](f)(P[])*(y-P[2])+f(P[]), x=-0.3..2.3, y=-0.3..2.3, style=surface, color="LightBlue"):
Point:=plottools:-sphere(P, 0.07, style=surface, color=red):
NormalLine:=plots:-spacecurve([P[1]+D[1](f)(P[])*t,P[2]+D[2](f)(P[])*t,P[3]-t], t=-1.7..1.7, color="Green", thickness=2):
plots:-display(Surface, Plane, Point, NormalLine, axes=normal);

If we express  x2(t)  from the first equation and substitute it into the second equation, we obtain a fourth-order linear equation with constant coefficients for the function  x1(t)  with a known right-hand side, that can be solved by known methods:

Eq1:=diff(x1(t),t,t)+b1*x1(t)-b2*x2(t)=f1(t):
Eq2:=diff(x2(t),t,t)+b3*x2(t)-b4*x2(t)=f2(t):
eval(Eq2, x2(t)=solve(Eq1,x2(t)));
collect(%, [diff(x1(t),t$4), diff(x1(t), t, t), x1(t)]);
% - `+`(op([1, 4..6], %));

Your example in Maple:

restart;
a:=unapply(rsolve({a(n)=-3*a(n-1)+a(n-2), a(0)=7, a(1)=2}, a(n)), n);  # General formula
seq(expand(a(n)), n=0..10);  # The first 11 terms

Statistics:-Histogram([2, 6], discrete=true, thickness=4, tickmarks = [[2=2*h^2, 6=6*h^2], default], view=[0..7, 0..0.55], axes=normal, axesfont=[times,roman,14]);

restart;
a:=100:  
b:=10:
k:=unapply(rsolve({k(n)=p*k(n-1), k(1)=b}, k(n)), n);    
x:=unapply(rsolve({x(n)=p*x(n-1)-k(n), x(1)=a-b}, x(n)), n);
z:=unapply(sum(x(n), n=1..k), k);
normal(z(30));
Optimization:-Maximize(%, p=0..2);

Addition: visual check.

plot(z(30), p=0..1, -50..300);

It's easier to use  plot  command:

plot([ [-1,4], [0,3], [1,4] ]);

1. Maple poorly works with sums with an undefined number of terms. If the number of terms is specific, then there is no problem:

restart;
is(sum(exp(x*i), i=1..10)>0)  assuming x::real;
                                                                                   true

2. In the second example, you made the mistake yourself:

is(ln(a+b)-ln(a)=ln(1+b/a)) assuming a>0, b>0;
                                                                                   true

Usually  return  is used only in procedures. In order to find the positions of all columns that have duplicate elements, you can simply save these positions as values of some indexed variable, and at the end convert these values to a set to exclude repetitions.

Example:

restart;
Dati:=LinearAlgebra:-RandomMatrix(5, 10, generator=1..10);
k:=0:
for j from 1 to 10 do
for i from 1 to 5 do
h:=Dati[i,j];
for x from 1 to 5 do        
if Dati[x,j] = h then
if x <> i then    k:=k+1; L[k]:=j;
end if; end if;
end do: end do: end do:
L:=convert(L, set);

Addition. Of course, this problem can be more effectively solved, for example, so

seq(`if`(numelems(Dati[..,k])<>nops(convert(Dati[..,k], set)), k, NULL), k=1..10);
                                        1, 3, 4, 5, 6, 9, 10