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

Example:

interface(rtablesize=infinity):
A:=LinearAlgebra:-RandomMatrix(30, 10, generator = 1 .. 9);
seq(convert(A[..,k], list), k=1..10);

In fact, your polynomial   r^2 – r*s + s^2  is a positive definite quadratic form. This means that the inequality   r^2 – r*s + s^2>0  holds for any values of the variables  r  and  s  that are not simultaneously equal to 0. Sylvester's criterion (see https://en.wikipedia.org/wiki/Sylvester%27s_criterion) checks the quadratic form of any number of variables whether it to be sign-definite. I did not find the appropriate procedure in Maple and therefore wrote my own procedure with the name Sylvester:

Sylvester:=proc(Q::polynom, Var::list)
local Coeff, n, A;
uses LinearAlgebra;
Coeff:=proc(P, Var, t)  # Extracts the coefficient in front of the monomial  t of a multivariate polynomial  P  of variables  Var
local L,H,i,k:
L:=[coeffs(P, Var, 'h')]: H:=[h]: k:=0:
for i from 1 to nops(H) do
if H[i]=t then k:=L[i] fi:
od:
k;
end proc;
n:=nops(Var);
A:=Matrix(n, (i,j)->`if`(i=j, Coeff(Q,Var,Var[i]*Var[j]), Coeff(Q,Var,Var[i]*Var[j])/2));
if `and`(seq(Determinant(A[1..k,1..k])>0, k=1..n)) then return
`Positive definite` elif
 `and`(seq(Determinant(-A[1..k,1..k])>0, k=1..n)) then return
`Negative definite` elif
 `and`(seq(Determinant(A[1..k,1..k])>=0, k=1..n)) then return
`Positive semidefinite` elif
`and`(seq(Determinant(-A[1..k,1..k])>=0, k=1..n)) then return 
`Negative semidefinite` else
`Sign-indefinite` fi;
end proc:

Examples of use.

Your example:

Sylvester(t^2-t*s+s^2, [t,s]);
                                                       Positive definite

More complicated examples:

Sylvester(2*x1^2-8*x1*x2-4*x1*x3+9*x2^2+19*x3^2, [x1,x2,x3]);
                                                        Positive definite

Sylvester(x1*x2+x2*x3-x1^2-x2^2-x3^2,[x1,x2,x3]);
                                                        Negative definite

restart;
a:=1000:  
b:=100:    
x[1]:=a-b:  k[1]:=b:
for i from 2 to 30 do
k[i]:=2*k[i-1];  x[i]:=2*x[i-1]-k[i];
od:
z=add(x[i], i=1..30);

                                                     z = -2040109466700

Addition. Of course, this problem is better solved using  rsolve  command. This allows us to obtain explicit formulas for the sequences  x(n)  and  k(n) :

restart;
a:=1000:  
b:=100:
k:=unapply(rsolve({k(n)=2*k(n-1), k(1)=b}, k(n)), n);    
x:=unapply(rsolve({x(n)=2*x(n-1)-k(n), x(1)=a-b}, x(n)), n);
z:=unapply(sum(x(n), n=1..k), k);
z(30);
z(1000);

restart;
local D:     
d := s*x*(E*K*q-K*r+K*sigma[1]+r*x)*
     (1-x*(E*K*q-K*r+K*sigma[1]+r*x)
     /(K*sigma[2]*L))/(K*sigma[2])+sigma[1]*x
     -x*(E*K*q-K*r+K*sigma[1]+r*x)/K:
L := ldegree(d,x):
P:=d/x^L:    
x^L * (add(``(simplify(coeff(P,x^k),size))*x^k, k=[3,2,1])+``(simplify(tcoeff(P,x),size)));
([A,B,C,D]=~ListTools:-Reverse(expand~([coeffs(op(2,%), x)])))[ ];

The issue is that  gcd  command only works for two polynomials. Therefore, the simplest solution for your example is to apply it 2 times:

P:=[x^2-7*x+10, x^2, x^2+2*x+1];
gcd(P[1], gcd(P[2], P[3]));
                                                                 1

If you have many polynomials, then you can use a simple user procedure gcdp :

gcdp:=proc(P::list(polynom))
if nops(P)=2 then return gcd(op(P)) else
gcd(P[1], gcdp(P[2..-1]));
fi;
end proc:

Example of use for your example:

gcdp(P);
                                                      1