I think Maple simply cannot cope with the symbolic (that is, exact) calculation of your integral. Therefore, I advise you to calculate it numerically, replacing infinity for example with the number 10. For comparison, I replaced infinity with the number 20. The results are identical. The reason is that the integrand very quickly converges to 0 (see its plot):

restart;

g(x,y):=min(3+(x-y)^2/10+abs(x+y)/sqrt(2),-abs(x-y)+7/sqrt(2));

q:=0; 
h(x,y):=Heaviside(g(x,y)-q);

p := Int(Int(h(x, y)*exp((-x^2-y^2)*(1/2))/(2*Pi), y = -10 .. 10), x = -10 .. 10);
evalf(%);

plot3d(h(x, y)*exp((-x^2-y^2)*(1/2))/(2*Pi), x=-5..5,y=-5..5);

p1 := Int(Int(h(x, y)*exp((-x^2-y^2)*(1/2))/(2*Pi), y = -20 .. 20), x = -20 .. 20);
evalf(%);

The number of subintervals is determined by the  partition  option. By default it is [5, 5] , that is, the interval for each variable is divided into 5 parts. Here are 2 examples with different values:

Download int.mw

When there is no idea how to symbolically represent a numeric expression, the following trick is useful:

identify(evalf((1/2)!));

You can easily do this using just 2 commands (in the title):

A:=DETools:-DEplot(diff(y(t),t)=(3-y(t))*(y(t)+1), y(t), t=0..5, y=-2..4, color=yellow, linecolor=white, [y(0)=4]):
B:=Student:-NumericalAnalysis:-Euler(diff(y(t),t)=(3-y(t))*(y(t)+1), y(0) = 4, t = 5, plotoptions=[color=red, symbolsize=12], output=plot):
plots:-display(A,B,scaling=constrained, size=[400,500]);

The nodal points in the Euler method are highlighted with small red circles (Maple does this automatically).

I can not confirm this. Everything works properly. Take a look:

restart;
F := x^2+y^2 = 1;
G := lhs(F) <= rhs(F);
is(eval(G, [x=0.5,y=0.6]));

                                       F := x^2+y^2 = 1
                                      G := x^2+y^2 <= 1
                                              true

Please submit the full text of the code in which you encountered this error "Error, cannot determine if this expression is true or false: x^2+y^2 <= 1"

Two syntax errors:

1. Should be  Pi  not  pi (it's just a symbol).

2. Should be  e(1)  instead of  e(T=1) .

Download 11.mw

Here is a simple procedure that implements your task. The default upper limit for the number of units to be checked is 5,000. You can easily change it in the procedure code if you wish.

Wrapped_prime:=proc(p::prime, N::posint:=5000)
local n, k, m0, m;
n:=length(p);
for k from 1 to N do
m0:=add(10^i, i=0..k-1);
m:=m0+10^k*p+10^(k+n)*m0;
if isprime(m) then return k fi;
od;
end proc:

Examples of use:

Wrapped_prime(59);
                                                  42

Wrapped_prime(7);
                                                   3

Wrapped_prime(2);
                                                  NULL

Wrapped_prime(71);
                                                   73

Wrapped_prime(167);
                                                  192
 

Here is another significantly more efficient way to generate primitive Pythagorean triples. This method is based on building a tree of such triples. At each step, the number of triples tripled. The simplest triple  <3, 4 ,5>  is the root of this tree. See this article  https://en.wikipedia.org/wiki/Tree_of_primitive_Pythagorean_triples  for details.

Below the code the procedure  PPT  (Primitive Pythagorean Triples)  and examples of use:

Download PPT.mw

Here is the simplest sorting algorithm called bubble sorting. See wiki for details.

restart;
BubbleSort:=proc(L::list(numeric))
local n, L1, k, m;
n:=nops(L);
L1:=L;
for k from 0 to n-2 do
for m from 1 to n-k-1 do
if L1[m]>L1[m+1] then L1:=subsop(m=L1[m+1],m+1=L1[m],L1) fi;
od; od;
L1;
end proc:

Example of use:

r:=rand(1..100):
L:=[seq(r(), i=1..20)];
BubbleSort(L);     

             L := [35, 15, 26, 100, 24, 8, 63, 78, 23, 73, 22, 32, 98, 9, 53, 3, 98, 69, 3, 73]
                   [3, 3, 8, 9, 15, 22, 23, 24, 26, 32, 35, 53, 63, 69, 73, 73, 78, 98, 98, 100]
 

Download plottest_new.mw

You can plot lines and points separately, and then all together using  plots:-display  command.

Example:

A:=plot(sin(x), x=0..2*Pi, color="DeepPink"):
B:=plot(cos(x), x=0..2*Pi, linestyle=3, color=blue):
A1:=plot([seq([x,sin(x)], x=0..6.29,evalf(Pi/2))], style=point, color="DeepPink", symbol=solidcircle, symbolsize=20):
B1:=plot([seq([x,cos(x)], x=0..6.29,evalf(Pi/3))], style=point, color=blue, symbol=solidbox, symbolsize=20):
plots:-display(A,B,A1,B1, size=[800,300], scaling=constrained);

To plot this contour, you do not need to solve any equation. You can simply substitute the value of this variable  Delta=-12.71  in the equation  f :

restart:
with(plots):

Omega:=0:lambda:=1:
gamma1:=8*Pi:
gamma2:=0.002*Pi:
x1:=100*Pi:
omega2:=200*Pi:
gamma0:=0.2*Pi:
G:=20*Pi:

lambda1:=(1/(2*Pi))*(G*Omega*gamma0/(2*(0.25*gamma0^2+Delta^2))):
lambda2:=(1/(2*Pi))*(-G*Omega*Delta/((0.25*gamma0^2+Delta^2))):
lambda3:=(1/(2*Pi))*gamma1+lambda1:
lambda4:=(1/(2*Pi))*(lambda*Delta-G^2*Delta/(0.25*gamma0^2+Delta^2)):
lambda5:=(1/(2*Pi))*(2*x1^2*omega2/((omega2^2+gamma2^2))):
f:=epsilon=(lambda1+sqrt(-lambda2^2+Y^2*(lambda3^2+(lambda4-lambda5*Y^2)^2)));

implicitplot3d(f, epsilon=0..20,Y = 0 .. 1,Delta=-40*Pi..40*Pi, labels=[epsilon,X,Delta],tickmarks=[3,3,3],style=surface);

implicitplot(eval(f,Delta=-12.71), epsilon=-20..20,Y = 0 .. 1, labels=[epsilon,Y], gridrefine=5);

L:=[.358700275060090779, [p[0] = .192413186080606, p[1] = 0.906594292940704e-1, p[2] = 0.677108912885780e-1, p[3] = 0.609551830556988e-1, p[4] = 0.589744573790909e-1, p[5] = 0.585737058072817e-1, p[6] = 0.589744573787748e-1, p[7] = 0.609551830550955e-1, p[8] = 0.677108912877626e-1, p[9] = 0.906594292931833e-1, p[10] = .192413186079858]]:
Data:=map(rhs, L[2]);
Statistics:-LineChart(Data);

It is better to solve your problems by writing mu and sigma down as procedures:

mu:=proc()
local k;
add(args[k],k=1..nargs)/nargs;
end proc:

sigma:=proc()
local k;
sqrt(add((args[k]-mu(args))^2, k=1..nargs)/(nargs-1));
end proc:

Examples of use:

mu($1..20);
sigma($1..20);

Codes of these procedures can be written more compactly if we use arrow notation and element-wise operator ~ :

mu:=()->`+`(args)/nargs:

sigma:=()->sqrt(add(([args]-~mu(args))^~2)/(nargs-1)):

Edit.

I think for a beginner the following code will be clearer. It does the same thing as  CartProd1  (of course for 3 lists only):

A:= [1,5,7]:
B:= [23,56,12]:
C:= [1,3]:
[seq(seq(seq([a,b,c],a=A),b=B),c=C)];